Data_analysis_amphibian_habitat_use
Julie Guenat
09/02/2023
R packages
library(tidyverse)
library(wesanderson)
library(pscl)
library(Matrix)
library(lme4)
library(DHARMa)
library(bbmle)
library(stats)
library(hrbrthemes)
library(ggpubr)
library(patchwork)
library(FactoMineR)
library(factoextra)
library(ape)
library(vegan)
library(corrplot)
library(ggordiplots)
library(jtools)
library(colorspace)
library(cowplot)
library(ggrepel)
library(ggsci)
library(ggtext)
library(ggthemes)
library(grid)
library(gridExtra)
library(rcartocolor)
library(scico)
library(showtext)
library(sf)
library(spatialreg)
library(spdep)
library(ade4)
library(adespatial)
library(adegraphics)
library(spdep)
library(sp)
library(zetadiv)
library(geosphere)
library(sjPlot)Raw reads data vizualisation
Relative reads abundance per Class
First, I intended to visualize the number of raw reads obtained for each Class retrieved using BATR01 primers. To do so, I loaded the tabfile I obtained after processing fasta.files with OBITools processing steps and removing the 1st PCR replicate that was not included for sequencing.
batr_raw<- read.table("C:/Users/jguenat/Documents/Paper_Master/2022_dataRaw_final/5_Raw_batr_WORep1.txt", header=T, sep="")
batr_raw<-as_tibble(batr_raw)I needed then to compute the real number of reads obtained in the different PCR replicated without the first one, since sequences were attributed to this replicate due to sequence leaking during sequencing process.
#total read count PCR replicate 1 included
Total_count<-sum(batr_raw$count)
Total_count## [1] 136704216#Total read count without PCR replicate 1
Total_count<-sum(batr_raw[,20:1075])
Total_count## [1] 131883594#Create a new column with counts corresponding to the total count without rep1
batr_raw<- batr_raw %>% mutate(count_real=rowSums(batr_raw[,20:1075]))The taxonomic levels at which sequences were assigned using the database, were: Order, family, genus and species. Since the taxonomic level of my interest is the class, I first needed to regroup the data per order and count the number of reads per order.
#determine the number of sequences per order
count_order<- batr_raw %>% group_by(order_name)%>%
summarise(sum(count_real))
#I exported this output as a file to be able to modify it.
##write.table(count_order, "C:/Users/jguenat/Documents/Paper_Master/2022_dataRaw_final/pieChart/1_Raw_Batr_count_order.txt", row.names = F)I attributed Class to each order using Excel by adding a column named “class”. I then loaded this file and computed the total read count for each class.
batr_count_class<-read.table("C:/Users/jguenat/Documents/Paper_Master/2022_dataRaw_final/pieChart/2_Raw_Batr_count_class.txt", header=T, sep="\t")
#count the number of sequences per Class and compute the sequences that it represents:
df<-batr_count_class %>% group_by(class)%>%
summarise(sum(sum.count_real.))%>%
mutate(rra=(`sum(sum.count_real.)`/sum(`sum(sum.count_real.)`))) %>%
mutate(rra_fig=c("0.27","0.38","0.02","0.27","0.06","<0.02"))Barplot of RRA per class
With ggplot2
cbPalette <- c("grey54", "indianred4",rep("grey54",4))
barplot_class <- ggplot(df, aes(x=rra, y=reorder(factor(class, labels=c("Actinopterygii", "Amphibia", "Aves", "Homo sapiens", "Unidentified","Other Mammalia")), rra), fill=class)) +
geom_bar(stat="identity", width = 0.9) +
labs(x="RRA") +
geom_text(aes(label = rra_fig), vjust = 0.45, hjust=-0.30) + #vjust = vertical; hjust = horizontal
scale_fill_manual(values = cbPalette, name="")+
theme_ipsum(grid="X") +
xlim(c(0,0.5)) +
theme(axis.text.x = element_text(size= 11, color= "black", vjust = -1),
axis.text.y = element_text(size= 11, color= "black"),
axis.title.x = element_text(size=12, vjust = -3, hjust = 0.5),
axis.title.y = element_blank(),
legend.position = "none")
barplot_class
Relative reads abundance per focal species
Then, I needed to extract the real number of reads for the species of interest and compute their relative reads abundance
sp<-batr_raw %>% filter(best_identity.Batr_3m_final>0.98 & species_name %in% c("Bufo bufo", "Hyla arborea","Lissotriton vulgaris", "Rana temporaria", "Pelophylax ridibundus", "Pelophylax bergeri")) %>% mutate(rra=count_real/sum(count_real))Since this study was conducted in two different reserves and that we estimated the densities of each of the focal amphibian species per reserve, I was interested in computing the relative read abundance per species and per reserve.
#Yverdon Reserve
yv <- sp %>% select(., contains(c("species_name", "pt1"))) %>%
mutate(numreads_yv = rowSums(.[,c(2:342)])) %>%
mutate(rra = numreads_yv / sum(numreads_yv)) %>%
mutate(reserve = "yv")#Gletterens Reserve
gl <- sp %>% select(., contains(c("species_name","pt2"))) %>%
mutate(numreads_gl = rowSums(.[,c(2:331)])) %>%
mutate(rra = numreads_gl / sum(numreads_gl))%>%
mutate(reserve="gle")Barplot of RRA per species and reserve
I plotted the relative read abundance of focal species per reserve. First, I needed to combine the rra from both reserve:
rra_sp<-bind_rows(gl, yv) %>%
select(-contains("pt"))
#View(rra_sp)Then, do the barplot per species and reserves:
palette2<-c("#9E6969", "#BD974C", "#BFC66E", "#91B3B8", "#7373B5", "#563B53")
barplot_sp<-ggplot(rra_sp, aes(y=rra, x= factor(reserve, labels = c("Gletterens", "Yverdon")), fill = factor(species_name, levels = c("Bufo bufo", "Hyla arborea", "Lissotriton vulgaris","Pelophylax bergeri", "Pelophylax ridibundus", "Rana temporaria")))) +
scale_fill_manual(values = palette2, "Species") +
geom_bar(stat = "identity", color = "black", width = 0.4)+
ylim(c(0,1)) + ylab("RRA")+ xlab("")+
geom_text(aes(label = round(rra,3)), hjust = -1)+
theme_ipsum(grid="X") +
theme(text= element_text(family ="A"),
axis.text.x= element_text(size = 11, color ="black"),
axis.text.y = element_text(size=11, color="black"),
axis.title.y = element_text(size=12, hjust= 0.5, vjust = 5),
axis.title.x = element_blank(),
legend.position = "none")
barplot_sp
Relation between RRA and estimated densities
Then, we wanted to investigate the relationship between the densities of amphibians and the rra.
Since we cannot differentiate Pelophylax ridibundus from bergeri morphologically speaking and we need to compute the densities of individuals retrieved with the prenuptial migration survey, we pondered the number of Pelophylax sp. individuals by the rra.
#yv ridi: rra= 0.048, berg: rra= 0.111, pelophylax sp. = 1
nridiyv<- 1*(0.048/(0.048+0.111))
nridiyv## [1] 0.3018868nbergyv<-1*(0.111/(0.111+0.048))
nbergyv## [1] 0.6981132#gle ridi: rra= 0.031, berg.= 0.139, pelophylax sp. = 2
nridigle<-2*(0.031/(0.031+0.139))
nridigle## [1] 0.3647059nberggle<-2*(0.139/(0.031+0.139))
nberggle## [1] 1.635294density<-read.table("C:/Users/jguenat/Documents/Paper_Master/Data_final/densities_sp_barriers_bergeri_pondere.txt", header=T, dec=",", sep="\t")To do so the first thing to do is to put rra and densities in a unique file
dens_rra<- left_join(rra_sp, density, by = c("species_name", "reserve"))Then I performed a linear regression on these term to test the statistical significance of the relationship between the density and the RRA.
mi<-lm(dens_rra$density~dens_rra$rra)
summary(mi)##
## Call:
## lm(formula = dens_rra$density ~ dens_rra$rra)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.7378 -4.6080 -3.6512 0.1592 29.8268
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.900 3.759 1.038 0.324
## dens_rra$rra 11.638 13.727 0.848 0.416
##
## Residual standard error: 10.33 on 10 degrees of freedom
## Multiple R-squared: 0.06706, Adjusted R-squared: -0.02623
## F-statistic: 0.7188 on 1 and 10 DF, p-value: 0.4164densi_plot<-ggplot(dens_rra, aes(x=density, y=rra))+
geom_point(aes(color=(species_name), shape=factor(reserve, labels = c("Yverdon", "Gletterens"))), size =4, alpha=1)+
scale_color_manual(values = palette2, "Species")+
scale_shape_manual(values= c(19,18), "Location")+
stat_smooth(method="lm", se=F, col="black") +
ylim(c(0:1))+ xlim(0,40)+
theme_minimal()+ ylab("RRA")+ xlab("Density [n/ha]")+
theme(axis.text = element_text(size=11, color = "black"),
axis.title.x = element_text(size=12),
axis.title.y = element_text(size = 12, vjust = -20),
axis.ticks = element_blank(),
legend.text = element_text(size=11, face="italic"),
legend.title = element_text(size=12),
legend.position = "none")+
stat_regline_equation(label.y = 1, label.x =0.01, aes(label = ..eq.label..)) +
stat_regline_equation(label.y = 0.90, label.x =0.01, aes(label = ..rr.label..)) +
guides(color = guide_legend(order = 1, keyheight = 1.5, keywidth = 1),
shape = guide_legend(order = 2, keyheight = 1.5, keywidth = 1))
densi_plot## `geom_smooth()` using formula = 'y ~ x'
It is also interesting to assess if the RRA is correlated to the number of sampling point where the species were found.
numpoint_plot<-ggplot(dens_rra, aes(x=pres_point, y=rra))+
geom_point(aes(color=(species_name), shape=factor(reserve, labels = c("Yverdon", "Gletterens"))), size =4, alpha=1)+
scale_color_manual(values = palette2, "Species")+
scale_shape_manual(values= c(19,18), "Location")+
stat_smooth(method="lm", se=F, col="black") +
ylim(c(0:1))+ xlim(0,25)+
theme_minimal()+ ylab("RRA")+ xlab("Number of detection points")+
theme(axis.text = element_text(size=11, color = "black"),
axis.title.x = element_text(size=12),
axis.title.y = element_text(size = 12, vjust = 5),
axis.ticks = element_blank(),
legend.text = element_text(size=11, face="italic"),
legend.title = element_text(size=12),
legend.position = "right")+
stat_regline_equation(label.y = 1, label.x =0.01, aes(label = ..eq.label..)) +
stat_regline_equation(label.y = 0.90, label.x =0.01, aes(label = ..rr.label..)) +
guides(color = guide_legend(order = 1, keyheight = 1.5, keywidth = 1),
shape = guide_legend(order = 2, keyheight = 1.5, keywidth = 1))
numpoint_plot## `geom_smooth()` using formula = 'y ~ x'
Then we statistically determine the relationship between the number of point where the spe. was detected and its rra.
ptdetec<-lm(dens_rra$pres_point~dens_rra$rra)
summary(ptdetec)##
## Call:
## lm(formula = dens_rra$pres_point ~ dens_rra$rra)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.8919 -3.9481 -0.6081 2.1735 13.0117
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.890 1.912 3.604 0.00482 **
## dens_rra$rra 1.661 6.982 0.238 0.81678
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.255 on 10 degrees of freedom
## Multiple R-squared: 0.005627, Adjusted R-squared: -0.09381
## F-statistic: 0.05659 on 1 and 10 DF, p-value: 0.8168Figure 2
To produce the fig. 2 of the manuscript, I combined the four plots above.
fig2<-(barplot_class + barplot_sp)/(densi_plot + numpoint_plot) + plot_layout(guides ='collect')# package patchwork
fig2_final<-fig2 + plot_annotation(tag_levels = 'A') &
theme(plot.tag = element_text(size = 12),
legend.spacing = unit(.5, "cm"),
legend.spacing.x = unit(0.5, "cm"))
fig2_final## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
Surface covered by the eDNA sampling
Since there is no relationship between RRA and estimated densities of the local amphibian species, I wanted to point out that the sampled surface using eDNA was largely smaller than the total surface of each reserve. Also the eDNA survey was not meant to monitor the density of each species, but rather investigate where these species occur.
So I wanted to compute the surface covered by the eDNA survey: A sampling point represents a circle of 5m diameter. The surface of a point is equal to pi* (rayon)^2
surface <- pi*(2.5)^2
round(surface, 2)## [1] 19.63We have 25 sampling point per reserve:
19.63*25# m^2## [1] 490.75Environmental data selection
We recorded 13 environmental variable to caracterize the amphibian habitat. However, we only have 50 sampling points and thus for further analyses we need to select the five most relevant variables. To do so, I performed a correlation plot to have an overview of the redunduncy of the variable and then selected among the least correlated variable, the most biologically relevant.
env<-read.table("C:/Users/jguenat/Documents/Paper_Master/2022_dataRaw_final/Env_data/final_env_Dec20.txt", header=T)
env<-as_tibble(env)
env.corr<-cor(env[,c(5:7, 9:18)])
corrplot(env.corr)
PCA(env[,c(5:7, 9:18)])
## **Results for the Principal Component Analysis (PCA)**
## The analysis was performed on 50 individuals, described by 13 variables
## *The results are available in the following objects:
##
## name description
## 1 "$eig" "eigenvalues"
## 2 "$var" "results for the variables"
## 3 "$var$coord" "coord. for the variables"
## 4 "$var$cor" "correlations variables - dimensions"
## 5 "$var$cos2" "cos2 for the variables"
## 6 "$var$contrib" "contributions of the variables"
## 7 "$ind" "results for the individuals"
## 8 "$ind$coord" "coord. for the individuals"
## 9 "$ind$cos2" "cos2 for the individuals"
## 10 "$ind$contrib" "contributions of the individuals"
## 11 "$call" "summary statistics"
## 12 "$call$centre" "mean of the variables"
## 13 "$call$ecart.type" "standard error of the variables"
## 14 "$call$row.w" "weights for the individuals"
## 15 "$call$col.w" "weights for the variables"Since there is a lot of variable, I decided to remove the proportion of each vegetation type.
env1.corr<-cor(env[, c(5:7, 9:11, 16:18)])
corrplot(env1.corr)
Then, I’ll check the continuous variables:
env1<-env[, c(5:7, 9:11, 16:18)]
as.tibble(env1)## # A tibble: 50 × 9
## X1_submerg X1_emerg X1_emerg_Land mean.water mean.mud forest.dist MeanTemp
## <int> <int> <int> <dbl> <dbl> <dbl> <dbl>
## 1 4 97 70 2.52 6.17 66.9 16.3
## 2 3 65 15 9.02 4.77 81.6 17.9
## 3 0 95 75 0.856 2.5 32.9 16.7
## 4 0 80 60 2.89 0.838 60.3 17.7
## 5 0 95 55 0.5 1.6 145. 16.8
## 6 1 70 60 2.67 3.37 95.3 19.3
## 7 1 93 35 4.66 1.89 58.1 16.6
## 8 1 60 45 0.856 1.7 140. 20.2
## 9 1 90 10 6.34 4.29 69.0 16.5
## 10 1 75 40 2.74 1.49 60.5 18.6
## # ℹ 40 more rows
## # ℹ 2 more variables: MinTemp <dbl>, MaxTemp <dbl>continuous<-select_if(env1, is.numeric)
#summary(continuous)# not same scale --> rescaling
#standardize the numeric variables:
env_scales<-env1 %>% mutate_if(is.numeric, funs(as.numeric(scale(.))))PCA(env_scales)
## **Results for the Principal Component Analysis (PCA)**
## The analysis was performed on 50 individuals, described by 9 variables
## *The results are available in the following objects:
##
## name description
## 1 "$eig" "eigenvalues"
## 2 "$var" "results for the variables"
## 3 "$var$coord" "coord. for the variables"
## 4 "$var$cor" "correlations variables - dimensions"
## 5 "$var$cos2" "cos2 for the variables"
## 6 "$var$contrib" "contributions of the variables"
## 7 "$ind" "results for the individuals"
## 8 "$ind$coord" "coord. for the individuals"
## 9 "$ind$cos2" "cos2 for the individuals"
## 10 "$ind$contrib" "contributions of the individuals"
## 11 "$call" "summary statistics"
## 12 "$call$centre" "mean of the variables"
## 13 "$call$ecart.type" "standard error of the variables"
## 14 "$call$row.w" "weights for the individuals"
## 15 "$call$col.w" "weights for the variables"par(mfrow=c(3,3))
hist(env_scales$X1_submerg, breaks = 30)
hist(env_scales$X1_emerg, breaks = 30)
hist(env_scales$X1_emerg_Land, breaks = 30)
hist(env_scales$mean.water, breaks = 30)
hist(env_scales$mean.mud, breaks = 30)
hist(env_scales$forest.dist, breaks = 30)
hist(env_scales$MeanTemp, breaks = 30)
hist(env_scales$MinTemp, breaks = 30)
hist(env_scales$MaxTemp, breaks = 30)
par(mfrow=c(3,3))
plot(env_scales$X1_submerg, breaks = 30)
plot(env_scales$X1_emerg, breaks = 30)
plot(env_scales$X1_emerg_Land, breaks = 30)
plot(env_scales$mean.water, breaks = 30)
plot(env_scales$mean.mud, breaks = 30)
plot(env_scales$forest.dist, breaks = 30)
plot(env_scales$MeanTemp, breaks = 30)
plot(env_scales$MinTemp, breaks = 30)
plot(env_scales$MaxTemp, breaks = 30)
pairs(env_scales)
From this I decided to reject Mean water depth, since it is quite correlated to the vegetation as well as the temperature. Moreover it is not necessarily biologically interesting. Indeed, water depth will determine the temperature but also the predation risk (as the more depths the more fishes) but I don’t have any other data to support this, so it is not interesting in my point of view.
I’ll keep the distance to the wintering habitats (estimation of the dispersion ability of sp.) and the mud depth (since important for the wintering habitat of Pelophylax sp. + spots for hunt or hind of sp.).
To determine which other environmental variables I’ll use, I did a last corr plot:
env2.corr<-cor(env_scales)
corrplot(env2.corr)
PCA(env_scales)
## **Results for the Principal Component Analysis (PCA)**
## The analysis was performed on 50 individuals, described by 9 variables
## *The results are available in the following objects:
##
## name description
## 1 "$eig" "eigenvalues"
## 2 "$var" "results for the variables"
## 3 "$var$coord" "coord. for the variables"
## 4 "$var$cor" "correlations variables - dimensions"
## 5 "$var$cos2" "cos2 for the variables"
## 6 "$var$contrib" "contributions of the variables"
## 7 "$ind" "results for the individuals"
## 8 "$ind$coord" "coord. for the individuals"
## 9 "$ind$cos2" "cos2 for the individuals"
## 10 "$ind$contrib" "contributions of the individuals"
## 11 "$call" "summary statistics"
## 12 "$call$centre" "mean of the variables"
## 13 "$call$ecart.type" "standard error of the variables"
## 14 "$call$row.w" "weights for the individuals"
## 15 "$call$col.w" "weights for the variables"I decided to keep: (1) Distance to the wintering habitat; (2) mud depth; (3) mean Water T°C; (4) proportion of emerged vegetation; (5) proportion of submerged vegetation; (6) proportion of emerged land.
env.reduced<-env[, c(1:7, 10:11, 16)]
##write.table(env.reduced,"C:/Users/jguenat/Documents/Paper_Master/2022_dataRaw_final/Env_data/2022final_envreduced.txt" ,row.names = F)
env.reduced$points<-as.character(env.reduced$points)
env.scaled<-env.reduced %>% mutate_if(is.numeric, funs(as.numeric(scale(.))))
##write.table(env.scaled,"C:/Users/jguenat/Documents/Paper_Master/2022_dataRaw_final/Env_data/2022final_envreduced_scaled.txt" ,row.names = F)Final Dataset
In order to analyse the effect of environmental factors on the probability of presence of amphibians sp., we first needed to combine the dataset with Pres_Abs of each species with the environmental variables dataset.
env<-read.table("C:/Users/jguenat/Documents/Paper_Master/2022_dataRaw_final/Env_data/2022final_envreduced.txt", header=T)
env<-as_tibble(env)
#scaled:
env2<-read.table("C:/Users/jguenat/Documents/Paper_Master/2022_dataRaw_final/Env_data/2022final_envreduced_scaled.txt", header=T)
env2<-as_tibble(env2)
batr<-read.table("C:/Users/jguenat/Documents/Paper_Master/2022_dataRaw_final/8_Batr_presabs.txt", header=T,sep="\t", check.names = F)
batr<-as_tibble(batr)
batr1<-batr %>% select(species_name, 11:60) %>%
pivot_longer(2:51, names_to="points", values_to="presabs")
env$points<-as.character(env$points)
evabund<-left_join(env, batr1, by="points")
#scaled:
env2$points<-as.character(env2$points)
evabund2<-left_join(env2, batr1, by="points")
#For convenience for further analyses, I will put each species as a column instead of in line:
dat1<-evabund2
bbufo<- dat1[dat1$species_name=="Bufo bufo",]
b<-bbufo %>% select(points, presabs) %>%
rename(Bbuf_pres=presabs)
dat1<-right_join(dat1, b)
Harbo<- dat1[dat1$species_name=="Hyla arborea",]
h<-Harbo %>% select(points, presabs) %>%
rename(Harbo_pres=presabs)
dat1<-right_join(dat1, h)
Lvulga<- dat1[dat1$species_name=="Lissotriton vulgaris",]
L<-Lvulga %>% select(points, presabs) %>%
rename(Lvulga_pres=presabs)
dat1<-right_join(dat1, L)
Pridi<- dat1[dat1$species_name=="Pelophylax ridibundus",]
P<-Pridi %>% select(points, presabs) %>%
rename(Pridi_pres=presabs)
dat1<-right_join(dat1, P)
Pberg<- dat1[dat1$species_name=="Pelophylax bergeri",]
pb<-Pberg %>% select(points, presabs) %>%
rename(Pberg_pres=presabs)
dat1<-right_join(dat1, pb)
Rtempo<- dat1[dat1$species_name=="Rana temporaria",]
R<- Rtempo %>% select(points, presabs) %>%
rename(Rtemp_pres=presabs)
dat1<-right_join(dat1, R)
#View(dat1)
x<- dat1 %>% select(!c(species_name, presabs)) %>% distinct()
#View(x)
##write.table(x,"C:/Users/jguenat/Documents/Paper_Master/2022_dataRaw_final/2022_FinalEnv_Batr_scaled.txt", row.names = F)Statistical analyses
data<- read.table("C:/Users/jguenat/Documents/Paper_Master/2022_dataRaw_final/2022_FinalEnv_Batr_scaled.txt", header=T, sep="")Zeta diversity
“Zeta diversity is a concept in ecology used to understand the turnover of species across multiple sites or sampling units. Unlike traditional diversity measures that typically focus on individual sites (alpha diversity) or pairs of sites (beta diversity), zeta diversity generalizes the concept to any number of sites, providing a more holistic view of species turnover across landscapes or communities.
Definition: Zeta diversity (ζ-diversity) measures the number of shared species among a given number of sites. For example, ζ₁ represents the total number of species observed in any single site, ζ₂ the number of species shared by any two sites, ζ₃ the number shared by any three sites, and so on.
Purpose: It provides insights into the scaling of biodiversity, showing how species turnover changes with the number of sites considered. It can help in understanding the spatial distribution of biodiversity, the effects of habitat fragmentation, and species co-occurrence patterns.
Calculation:
ζ₁ is simply the average species richness across sites. ζ₂ is the average number of species shared between any pair of sites. ζ₃ is the average number of species shared among any three sites, and so on. The values can be calculated for k sites (ζₖ) and usually decrease as k increases, reflecting the likelihood that more sites share fewer species.
Zeta diversity is implemented in R using the zetadiv package.
Then, we need to prepare the data to compute the zeta diversity. For that we need one matrix with the species pres/abs and a matrix with the geographical distances.”
Yverdon
select the data:
data2<- read.table("C:/Users/jguenat/Documents/Paper_Master/2022_dataRaw_final/2024_coordadded.txt", header=T, sep="")
data_1 <- data2 %>% filter(Location == "Cheseaux")
xy = data_1[, c("x", "y")]
xy_coord <- st_as_sf(xy, coords = c("x", "y"), crs = 4326)
coordinates_utm <- st_transform(xy_coord, crs = 32632)
xy <- as.data.frame(st_coordinates(coordinates_utm))
occurences = data_1[, 11:16]for each order (2, 3, 5, 10), a linear regression and a gam are computed by default distance = Euclidean.
set.seed(1)
zeta.ddecay.lm.fine2 <- Zeta.ddecay(xy, occurences, sam = 1000, order = 2, confint.level = 0.95, method = "mean", plot = T)
set.seed(1)
zeta.ddecay.lm.fine3 <- Zeta.ddecay(xy, occurences, sam = 1000, order = 3, confint.level = 0.95, method = "mean", plot = T)
set.seed(1)
zeta.ddecay.lm.fine5 <- Zeta.ddecay(xy, occurences, sam = 1000, order = 5, confint.level = 0.95, method = "mean", plot = T)
set.seed(1)
zeta.ddecay.lm.fine10 <- Zeta.ddecay(xy, occurences, sam = 1000, order = 10, confint.level = 0.95, method = "mean", plot = T)
Zeta 2
To be able to modify the graphs, we need to extract the information:
distance_z2 <- zeta.ddecay.lm.fine2$distance
zeta2<- zeta.ddecay.lm.fine2$zeta.val
data_plot_z2_yv <- data.frame(distance = distance_z2, zeta2= zeta2)
coeff_z2_yv <- coef(zeta.ddecay.lm.fine2$reg)
intercept_z2_yv <- coeff_z2_yv[1]
slope_z2_yv <- coeff_z2_yv[2]
#Extract predicted values and standar errors:
predictions_z2_yv <- predict(zeta.ddecay.lm.fine2$reg, newdata = data_plot_z2_yv, se.fit = TRUE)
# Extract predicted values and standard errors
data_plot_z2_yv$predicted <- predictions_z2_yv$fit
data_plot_z2_yv$se <- predictions_z2_yv$se.fit
# Calculate confidence intervals
data_plot_z2_yv$lower <- data_plot_z2_yv$predicted - 1.96 * data_plot_z2_yv$se
data_plot_z2_yv$upper <- data_plot_z2_yv$predicted + 1.96 * data_plot_z2_yv$seGraphical representation:
zeta2_yv <- ggplot(data_plot_z2_yv, aes(x = distance, y = zeta2)) +
geom_point(color = "black", size = 2, alpha = 0.25) +
geom_line(aes(y = predicted), color = "black", size = 1) + # Predicted values line
geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.35, fill = "darkgrey") +
labs(x = "Distance between two sites (m)",
y = expression("Number of species shared by two sites (ζ"[2]*")"),
title = "Yverdon") +
theme_light() +
theme(axis.title.x = element_text(size = 11),
axis.title.y = element_text(size = 11),
plot.title = element_text(hjust = 0.5, size = 12))## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.Zeta 3
distance_z3 <- zeta.ddecay.lm.fine3$distance
zeta3<- zeta.ddecay.lm.fine3$zeta.val
data_plot_z3_yv <- data.frame(distance = distance_z3, zeta3= zeta3)
coeff_z3_yv <- coef(zeta.ddecay.lm.fine3$reg)
intercept_z3_yv <- coeff_z3_yv[1]
slope_z3_yv <- coeff_z3_yv[2]
#Extract predicted values and standar errors:
predictions_z3_yv <- predict(zeta.ddecay.lm.fine3$reg, newdata = data_plot_z3_yv, se.fit = TRUE)
# Extract predicted values and standard errors
data_plot_z3_yv$predicted <- predictions_z3_yv$fit
data_plot_z3_yv$se <- predictions_z3_yv$se.fit
# Calculate confidence intervals
data_plot_z3_yv$lower <- data_plot_z3_yv$predicted - 1.96 * data_plot_z3_yv$se
data_plot_z3_yv$upper <- data_plot_z3_yv$predicted + 1.96 * data_plot_z3_yv$se
zeta3_yv <- ggplot(data_plot_z3_yv, aes(x = distance, y = zeta3)) +
geom_point(color = "black", size = 2, alpha = 0.25) +
geom_line(aes(y = predicted), color = "black", size = 1) + # Predicted values line
geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.35, fill = "darkgrey") +
labs(x = "Mean distance between three sites (m)",
y = expression("Number of species shared by three sites (ζ"[3]*")"),
title = "Yverdon") +
theme_light() +
theme(axis.title.x = element_text(size = 11),
axis.title.y = element_text(size = 11),
plot.title = element_text(hjust = 0.5, size = 12))Zeta 5
distance_z5 <- zeta.ddecay.lm.fine5$distance
zeta5<- zeta.ddecay.lm.fine5$zeta.val
data_plot_z5_yv <- data.frame(distance = distance_z5, zeta5= zeta5)
coeff_z5_yv <- coef(zeta.ddecay.lm.fine5$reg)
intercept_z5_yv <- coeff_z5_yv[1]
slope_z5_yv <- coeff_z5_yv[2]
#Extract predicted values and standar errors:
predictions_z5_yv <- predict(zeta.ddecay.lm.fine5$reg, newdata = data_plot_z5_yv, se.fit = TRUE)
# Extract predicted values and standard errors
data_plot_z5_yv$predicted <- predictions_z5_yv$fit
data_plot_z5_yv$se <- predictions_z5_yv$se.fit
# Calculate confidence intervals
data_plot_z5_yv$lower <- data_plot_z5_yv$predicted - 1.96 * data_plot_z5_yv$se
data_plot_z5_yv$upper <- data_plot_z5_yv$predicted + 1.96 * data_plot_z5_yv$se
zeta5_yv <- ggplot(data_plot_z5_yv, aes(x = distance, y = zeta5)) +
geom_point(color = "black", size = 2, alpha = 0.25) +
geom_line(aes(y = predicted), color = "black", size = 1) + # Predicted values line
geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.35, fill = "darkgrey") +
labs(x = "Mean distance between five sites (m)",
y = expression("Number of species shared by five sites (ζ"[5]*")"),
title = "Yverdon") +
theme_light() +
theme(axis.title.x = element_text(size = 11),
axis.title.y = element_text(size = 11),
plot.title = element_text(hjust = 0.5, size = 12))Gletterens
data_2 <- data2 %>% filter(Location == "Glettrens")
xy = data_2[, c("x", "y")]
xy_coord <- st_as_sf(xy, coords = c("x", "y"), crs = 4326)
coordinates_utm <- st_transform(xy_coord, crs = 32632)
xy <- as.data.frame(st_coordinates(coordinates_utm))
occurences = data_2[, 11:16]for each order (2, 3, 5, 10), a linear regression and a gam are computed by default distance = Euclidean
set.seed(1)
zeta.ddecay.lm.fine2 <- Zeta.ddecay(xy, occurences, sam = 1000, order = 2, confint.level = 0.95, method = "mean", plot = T)
set.seed(1)
zeta.ddecay.lm.fine3 <- Zeta.ddecay(xy, occurences, sam = 1000, order = 3, confint.level = 0.95, method = "mean", plot = T)
set.seed(1)
zeta.ddecay.lm.fine5 <- Zeta.ddecay(xy, occurences, sam = 1000, order = 5, confint.level = 0.95, method = "mean", plot = T)
set.seed(1)
zeta.ddecay.lm.fine10 <- Zeta.ddecay(xy, occurences, sam = 1000, order = 10, confint.level = 0.95, method = "mean", plot = T)
Zeta 2
distance_z2 <- zeta.ddecay.lm.fine2$distance
zeta2<- zeta.ddecay.lm.fine2$zeta.val
data_plot_z2_gle <- data.frame(distance = distance_z2, zeta2= zeta2)
coeff_z2_gle <- coef(zeta.ddecay.lm.fine2$reg)
intercept_z2_gle <- coeff_z2_gle[1]
slope_z2_gle <- coeff_z2_gle[2]
#Extract predicted values and standar errors:
predictions_z2_gle <- predict(zeta.ddecay.lm.fine2$reg, newdata = data_plot_z2_gle, se.fit = TRUE)
# Extract predicted values and standard errors
data_plot_z2_gle$predicted <- predictions_z2_gle$fit
data_plot_z2_gle$se <- predictions_z2_gle$se.fit
# Calculate confidence intervals
data_plot_z2_gle$lower <- data_plot_z2_gle$predicted - 1.96 * data_plot_z2_gle$se
data_plot_z2_gle$upper <- data_plot_z2_gle$predicted + 1.96 * data_plot_z2_gle$se
zeta2_gle <- ggplot(data_plot_z2_gle, aes(x = distance, y = zeta2)) +
geom_point(color = "black", size = 2, alpha = 0.25) +
geom_line(aes(y = predicted), color = "black", size = 1) + # Predicted values line
geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.35, fill = "darkgrey") +
labs(x = "Distance between two sites (m)",
y = expression("Number of species shared by two sites (ζ"[2]*")"),
title = "Gletterens") +
theme_light() +
theme(axis.title.x = element_text(size = 11),
axis.title.y = element_text(size = 11),
plot.title = element_text(hjust = 0.5, size = 12))Zeta 3
distance_z3 <- zeta.ddecay.lm.fine3$distance
zeta3<- zeta.ddecay.lm.fine3$zeta.val
data_plot_z3_gle <- data.frame(distance = distance_z3, zeta3= zeta3)
coeff_z3_gle <- coef(zeta.ddecay.lm.fine3$reg)
intercept_z3_gle <- coeff_z3_gle[1]
slope_z3_gle <- coeff_z3_gle[2]
#Extract predicted values and standar errors:
predictions_z3_gle <- predict(zeta.ddecay.lm.fine3$reg, newdata = data_plot_z3_gle, se.fit = TRUE)
# Extract predicted values and standard errors
data_plot_z3_gle$predicted <- predictions_z3_gle$fit
data_plot_z3_gle$se <- predictions_z3_gle$se.fit
# Calculate confidence intervals
data_plot_z3_gle$lower <- data_plot_z3_gle$predicted - 1.96 * data_plot_z3_gle$se
data_plot_z3_gle$upper <- data_plot_z3_gle$predicted + 1.96 * data_plot_z3_gle$se
zeta3_gle <- ggplot(data_plot_z3_gle, aes(x = distance, y = zeta3)) +
geom_point(color = "black", size = 2, alpha = 0.25) +
geom_line(aes(y = predicted), color = "black", size = 1) + # Predicted values line
geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.35, fill = "darkgrey") +
labs(x = "Mean distance between three sites (m)",
y = expression("Number of species shared by three sites (ζ"[3]*")"),
title = "Gleterrens") +
theme_light() +
theme(axis.title.x = element_text(size = 11),
axis.title.y = element_text(size = 11),
plot.title = element_text(hjust = 0.5, size = 12))Zeta 5
distance_z5 <- zeta.ddecay.lm.fine5$distance
zeta5<- zeta.ddecay.lm.fine5$zeta.val
data_plot_z5_gle <- data.frame(distance = distance_z5, zeta5= zeta5)
coeff_z5_gle <- coef(zeta.ddecay.lm.fine5$reg)
intercept_z5_gle <- coeff_z5_gle[1]
slope_z5_gle <- coeff_z5_gle[2]
#Extract predicted values and standar errors:
predictions_z5_gle <- predict(zeta.ddecay.lm.fine5$reg, newdata = data_plot_z5_gle, se.fit = TRUE)
# Extract predicted values and standard errors
data_plot_z5_gle$predicted <- predictions_z5_gle$fit
data_plot_z5_gle$se <- predictions_z5_gle$se.fit
# Calculate confidence intervals
data_plot_z5_gle$lower <- data_plot_z5_gle$predicted - 1.96 * data_plot_z5_gle$se
data_plot_z5_gle$upper <- data_plot_z5_gle$predicted + 1.96 * data_plot_z5_gle$se
zeta5_gle <- ggplot(data_plot_z5_gle, aes(x = distance, y = zeta5)) +
geom_point(color = "black", size = 2, alpha = 0.25) +
geom_line(aes(y = predicted), color = "black", size = 1) + # Predicted values line
geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.35, fill = "darkgrey") +
labs(x = "Mean distance between five sites (m)",
y = expression("Number of species shared by five sites (ζ"[5]*")"),
title = "Gletterens") +
theme_light() +
theme(axis.title.x = element_text(size = 11),
axis.title.y = element_text(size = 11),
plot.title = element_text(hjust = 0.5, size = 12))Plot zeta diversity
Fig_S4<-(zeta2_yv | zeta3_yv | zeta5_yv) / (zeta2_gle | zeta3_gle | zeta5_gle)
Fig_S4 + plot_annotation(tag_levels = 'A')
db-RDA replacing 0 by 0.001
We asked if the species assembly could be explained by the habitat characteristics proxied by the environmental variables we measured. We performed a distance-based redundancy analysis.
I wanted to have a first look at how the proportion of emerged, submerged vegetation and emerged land affect the species assembly, to determine if these variables are all important ones.
dat<-data
#dat[dat == 0]<- 0.001
sp=dat[,11:16]
sp001<- sp + 0.001
env=dat[,5:10]
#sp.dist001<-vegdist(sp001, method = "jaccard")
colnames(sp)<-c("Bb","Ha","Lv","Pr","Pb","Rt")
colnames(sp001)<-c("Bb","Ha","Lv","Pr","Pb","Rt")
colnames(env)<-c("SV", "EV", "EL", "MD", "FD", "WT")
rankindex(env, sp001, indices = c("jaccard", "bray"), stepacross= FALSE, method = "spearman")## jaccard bray
## -0.03328199 -0.03293306#Capscale function from vegan package
dbRDA = capscale(sp001 ~ SV + EV + EL + MD + FD + WT,
env, dist="jaccard")
summary(dbRDA)##
## Call:
## capscale(formula = sp001 ~ SV + EV + EL + MD + FD + WT, data = env, distance = "jaccard")
##
## Partitioning of squared Jaccard distance:
## Inertia Proportion
## Total 16.671 1.000
## Constrained 2.268 0.136
## Unconstrained 14.403 0.864
##
## Eigenvalues, and their contribution to the squared Jaccard distance
##
## Importance of components:
## CAP1 CAP2 CAP3 CAP4 CAP5 CAP6 MDS1
## Eigenvalue 0.99265 0.64674 0.35951 0.24218 0.023288 0.0031386 5.7253
## Proportion Explained 0.05955 0.03880 0.02157 0.01453 0.001397 0.0001883 0.3434
## Cumulative Proportion 0.05955 0.09834 0.11991 0.13443 0.135831 0.1360188 0.4795
## MDS2 MDS3 MDS4 MDS5 MDS6 MDS7 MDS8
## Eigenvalue 3.0656 2.2314 1.8295 0.71615 0.54941 0.2701 0.0155887
## Proportion Explained 0.1839 0.1339 0.1097 0.04296 0.03296 0.0162 0.0009351
## Cumulative Proportion 0.6634 0.7972 0.9069 0.94991 0.98286 0.9991 1.0000000
##
## Accumulated constrained eigenvalues
## Importance of components:
## CAP1 CAP2 CAP3 CAP4 CAP5 CAP6
## Eigenvalue 0.9926 0.6467 0.3595 0.2422 0.02329 0.003139
## Proportion Explained 0.4378 0.2852 0.1586 0.1068 0.01027 0.001384
## Cumulative Proportion 0.4378 0.7230 0.8815 0.9883 0.99862 1.000000RsquareAdj(dbRDA)## $r.squared
## [1] 0.1360188
##
## $adj.r.squared
## [1] 0.01546327Since the relation between the probability of presences of species and the environmental factors composing their habitat might not be linear (i.e. species might have an unimodal distribution, with an optimum and some tolerance to the variation), we tested the effect on the model of the addition of quadratic terms of all environmental variables. We used the polyvar() function proposed in “Numerical Ecology with R” from Bocard, Gillet & Legendre (2018).
# adding the quadratic effect of environmental variables:
env.square<- polyvars(env, degr = 2)
rankindex(env.square, sp001, indices = c("jaccard", "bray"), stepacross= FALSE, method = "spearman")## jaccard bray
## -0.04761048 -0.04725647dbRDA.square<- capscale(sp001 ~ ., env.square, distance = "jaccard")
summary(dbRDA.square)##
## Call:
## capscale(formula = sp001 ~ SV.1 + SV.2 + EV.1 + EV.2 + EL.1 + EL.2 + MD.1 + MD.2 + FD.1 + FD.2 + WT.1 + WT.2, data = env.square, distance = "jaccard")
##
## Partitioning of squared Jaccard distance:
## Inertia Proportion
## Total 16.671 1.000
## Constrained 4.568 0.274
## Unconstrained 12.102 0.726
##
## Eigenvalues, and their contribution to the squared Jaccard distance
##
## Importance of components:
## CAP1 CAP2 CAP3 CAP4 CAP5 CAP6
## Eigenvalue 2.1291 1.33975 0.52667 0.39283 0.123676 0.038943
## Proportion Explained 0.1277 0.08037 0.03159 0.02356 0.007419 0.002336
## Cumulative Proportion 0.1277 0.20808 0.23967 0.26324 0.270657 0.272993
## CAP7 CAP8 MDS1 MDS2 MDS3 MDS4 MDS5
## Eigenvalue 0.0155790 0.0018441 4.5776 2.8056 1.8136 1.45259 0.65768
## Proportion Explained 0.0009345 0.0001106 0.2746 0.1683 0.1088 0.08714 0.03945
## Cumulative Proportion 0.2739273 0.2740379 0.5486 0.7169 0.8257 0.91286 0.95231
## MDS6 MDS7 MDS8
## Eigenvalue 0.53390 0.2468 0.0143806
## Proportion Explained 0.03203 0.0148 0.0008626
## Cumulative Proportion 0.98433 0.9991 1.0000000
##
## Accumulated constrained eigenvalues
## Importance of components:
## CAP1 CAP2 CAP3 CAP4 CAP5 CAP6 CAP7
## Eigenvalue 2.129 1.3398 0.5267 0.39283 0.12368 0.038943 0.01558
## Proportion Explained 0.466 0.2933 0.1153 0.08599 0.02707 0.008524 0.00341
## Cumulative Proportion 0.466 0.7593 0.8746 0.96059 0.98766 0.996186 0.99960
## CAP8
## Eigenvalue 0.0018441
## Proportion Explained 0.0004037
## Cumulative Proportion 1.0000000RsquareAdj(dbRDA.square)## $r.squared
## [1] 0.2740379
##
## $adj.r.squared
## [1] 0.03859079The adjusted R.squared is way better with the quadratic effect of the environmental variables, thus we decided to keep them.
Ordination plot of db-RDA
First, the triplot for the db-RDA without quadratic effect:
shapes = c(22, 23, 24)
shapes <- shapes[as.factor(dat$environment)]
colors <- c("red", "darkgreen", "blue")
colors <- colors[as.factor(dat$environment)]
Names = paste("N", 1:4, sep=" = ")
plot(dbRDA, type= "n", scaling = 3, las=1)
points(dbRDA, display = "species", pch=20, cex = 1.5, col="#ff7f00", scaling=3)
text(dbRDA, display = "species", scaling= 3, labels = colnames(sp), cex=0.8, pos= 3)
points(dbRDA, display = "sites", pch= shapes, cex= 1.2, col= colors, scaling=3)
text(dbRDA, scaling= 3, display="bp", col="#0868ac", cex= 0.8)
"#0868ac"## [1] "#0868ac"# Legend
legend("topright", legend = levels(as.factor(dat$environment)),
col = c("red", "darkgreen", "blue"),
pch = c(22, 23, 24) )
Significance of db-RDA
We tested whether the influences of the environmental variables on the species assemblage were significant by performing a permanova (i.e., permutational multivariate analysis of variance) using the function adonis2 from the package ‘vegan’.
#PERMANOVA!!!!
adonis2(sp001 ~ SV + EV + EL + MD + FD + WT,
env, dist="jaccard", permutations = 999)## Permutation test for adonis under reduced model
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = sp001 ~ SV + EV + EL + MD + FD + WT, data = env, permutations = 999, dist = "jaccard")
## Df SumOfSqs R2 F Pr(>F)
## Model 6 1.9105 0.13564 1.1246 0.355
## Residual 43 12.1746 0.86436
## Total 49 14.0851 1.00000anova(dbRDA, permutations=how(nperm = 999))## Permutation test for capscale under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model: capscale(formula = sp001 ~ SV + EV + EL + MD + FD + WT, data = env, distance = "jaccard")
## Df SumOfSqs F Pr(>F)
## Model 6 2.2675 1.1283 0.295
## Residual 43 14.4030#with quadratic effect
anova(dbRDA.square, permutations=how(nperm = 999))## Permutation test for capscale under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model: capscale(formula = sp001 ~ SV.1 + SV.2 + EV.1 + EV.2 + EL.1 + EL.2 + MD.1 + MD.2 + FD.1 + FD.2 + WT.1 + WT.2, data = env.square, distance = "jaccard")
## Df SumOfSqs F Pr(>F)
## Model 12 4.5684 1.1639 0.226
## Residual 37 12.1022adonis2(sp001 ~ .,
env.square, dist="jaccard", permutations = 999)## Permutation test for adonis under reduced model
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = sp001 ~ ., data = env.square, permutations = 999, dist = "jaccard")
## Df SumOfSqs R2 F Pr(>F)
## Model 12 3.9623 0.28131 1.2069 0.236
## Residual 37 10.1228 0.71869
## Total 49 14.0851 1.00000##RDA
dat<-data
sp=dat[,11:16]
sp001<- sp + 0.001
env=dat[,5:10]
colnames(sp)<-c("Bb","Ha","Lv","Pr","Pb","Rt")
colnames(sp001)<-c("Bb","Ha","Lv","Pr","Pb","Rt")
colnames(env)<-c("SV", "EV", "EL", "MD", "FD", "WT")
#RDA function from vegan package
RDA = rda(sp001 ~ SV + EV + EL + MD + FD + WT,
env)
summary(RDA)##
## Call:
## rda(formula = sp001 ~ SV + EV + EL + MD + FD + WT, data = env)
##
## Partitioning of variance:
## Inertia Proportion
## Total 1.0898 1.0000
## Constrained 0.1471 0.1349
## Unconstrained 0.9427 0.8651
##
## Eigenvalues, and their contribution to the variance
##
## Importance of components:
## RDA1 RDA2 RDA3 RDA4 RDA5 RDA6
## Eigenvalue 0.06123 0.03067 0.02875 0.02468 0.0009816 0.0007378
## Proportion Explained 0.05619 0.02814 0.02639 0.02264 0.0009008 0.0006770
## Cumulative Proportion 0.05619 0.08433 0.11072 0.13336 0.1342601 0.1349371
## PC1 PC2 PC3 PC4 PC5 PC6
## Eigenvalue 0.3472 0.2523 0.1616 0.08730 0.05446 0.03985
## Proportion Explained 0.3186 0.2315 0.1483 0.08011 0.04997 0.03657
## Cumulative Proportion 0.4536 0.6850 0.8333 0.91346 0.96343 1.00000
##
## Accumulated constrained eigenvalues
## Importance of components:
## RDA1 RDA2 RDA3 RDA4 RDA5 RDA6
## Eigenvalue 0.06123 0.03067 0.02875 0.02468 0.0009816 0.0007378
## Proportion Explained 0.41640 0.20856 0.19554 0.16781 0.0066754 0.0050175
## Cumulative Proportion 0.41640 0.62496 0.82050 0.98831 0.9949825 1.0000000RsquareAdj(RDA)## $r.squared
## [1] 0.1349371
##
## $adj.r.squared
## [1] 0.01423067# adding the quadratic effect of environmental variables:
env.square<- polyvars(env, degr = 2)
RDA.square<- rda(sp001 ~ ., env.square)
summary(RDA.square)##
## Call:
## rda(formula = sp001 ~ SV.1 + SV.2 + EV.1 + EV.2 + EL.1 + EL.2 + MD.1 + MD.2 + FD.1 + FD.2 + WT.1 + WT.2, data = env.square)
##
## Partitioning of variance:
## Inertia Proportion
## Total 1.0898 1.000
## Constrained 0.2648 0.243
## Unconstrained 0.8250 0.757
##
## Eigenvalues, and their contribution to the variance
##
## Importance of components:
## RDA1 RDA2 RDA3 RDA4 RDA5 RDA6 PC1
## Eigenvalue 0.1119 0.06855 0.03992 0.02988 0.01117 0.003431 0.2891
## Proportion Explained 0.1027 0.06290 0.03663 0.02741 0.01025 0.003148 0.2653
## Cumulative Proportion 0.1027 0.16558 0.20221 0.22962 0.23987 0.243021 0.5083
## PC2 PC3 PC4 PC5 PC6
## Eigenvalue 0.2232 0.1507 0.07658 0.05212 0.03323
## Proportion Explained 0.2048 0.1383 0.07027 0.04782 0.03049
## Cumulative Proportion 0.7132 0.8514 0.92169 0.96951 1.00000
##
## Accumulated constrained eigenvalues
## Importance of components:
## RDA1 RDA2 RDA3 RDA4 RDA5 RDA6
## Eigenvalue 0.1119 0.06855 0.03992 0.02988 0.01117 0.003431
## Proportion Explained 0.4225 0.25884 0.15071 0.11280 0.04218 0.012956
## Cumulative Proportion 0.4225 0.68135 0.83206 0.94487 0.98704 1.000000RsquareAdj(RDA.square)## $r.squared
## [1] 0.2430214
##
## $adj.r.squared
## [1] -0.002485127# with squared effect:
Names = paste("N", 1:4, sep=" = ")
plot(dbRDA.square, type= "n", scaling = 2, las=1, xlab = "CAP 1 (12.77%)", ylab= "CAP 2 (8.04%)", ylim = c(-2.5, 2.5), xlim = c(-4, 4))
points(dbRDA.square, display = "species", pch=19, cex = 1.5, col=palette2, scaling=3)
text(dbRDA.square, display = "species", scaling= 3, labels = colnames(sp), cex=0.8, pos= 2)
points(dbRDA.square, display = "sites", pch = 20, cex= 1, col= "black", scaling=3)
text(dbRDA.square, scaling= 3, display="bp", col="#0868ac", cex= 0.8)
Generalized Linear Models
We aimed at characterizing the breeding habitat of each of the local amphibian species. To this end, we tested the effect of the recorded environmental variables on the probability of presence of each amphibian species. We performed separate binomial generalized linear models (GLMs) for each species, using the presence-absence of the focal species as response variable and the environmental variables as explanatory variables.
based on the graphs chunk 45, we can observe that the distribution of points does not have a linear pattern. Thus, I’ll add the quadratic effect of these variables in the following models.
Bufo bufo
length(data$points[data$Bbuf_pres==1])## [1] 19We wanted to analyze the effect of other amphibians species on the probability of presence of the focal sp. Since there is 5 other species and that we are already testing 6 environmental variable, we would end with 11 explanatory variable that is way too much for 50 sampling points. What we did then, is to first do a PCA on the presence of other amphibian sp and then use this as an explanatory variable in the model.
pcaB<-prcomp(data[,c(12:16)])
pcaB## Standard deviations (1, .., p=5):
## [1] 0.6281888 0.4309699 0.3815493 0.2621920 0.2338955
##
## Rotation (n x k) = (5 x 5):
## PC1 PC2 PC3 PC4 PC5
## Harbo_pres 0.1419059 0.10666751 -0.200815801 0.7733082 -0.57458873
## Lvulga_pres -0.1818155 -0.95494934 -0.006808097 0.2211551 0.07783914
## Pridi_pres -0.1483413 0.15594608 -0.672260323 0.3058557 0.63890062
## Pberg_pres -0.6728335 0.01282751 -0.453510094 -0.3459734 -0.47091568
## Rtemp_pres 0.6870923 -0.22849523 -0.549564661 -0.3739513 -0.18393825str(pcaB)## List of 5
## $ sdev : num [1:5] 0.628 0.431 0.382 0.262 0.234
## $ rotation: num [1:5, 1:5] 0.142 -0.182 -0.148 -0.673 0.687 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:5] "Harbo_pres" "Lvulga_pres" "Pridi_pres" "Pberg_pres" ...
## .. ..$ : chr [1:5] "PC1" "PC2" "PC3" "PC4" ...
## $ center : Named num [1:5] 0.08 0.24 0.12 0.34 0.56
## ..- attr(*, "names")= chr [1:5] "Harbo_pres" "Lvulga_pres" "Pridi_pres" "Pberg_pres" ...
## $ scale : logi FALSE
## $ x : num [1:50, 1:5] -0.927 -0.779 0.581 0.723 -0.106 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : NULL
## .. ..$ : chr [1:5] "PC1" "PC2" "PC3" "PC4" ...
## - attr(*, "class")= chr "prcomp"PCA(data[,c(12:16)])
## **Results for the Principal Component Analysis (PCA)**
## The analysis was performed on 50 individuals, described by 5 variables
## *The results are available in the following objects:
##
## name description
## 1 "$eig" "eigenvalues"
## 2 "$var" "results for the variables"
## 3 "$var$coord" "coord. for the variables"
## 4 "$var$cor" "correlations variables - dimensions"
## 5 "$var$cos2" "cos2 for the variables"
## 6 "$var$contrib" "contributions of the variables"
## 7 "$ind" "results for the individuals"
## 8 "$ind$coord" "coord. for the individuals"
## 9 "$ind$cos2" "cos2 for the individuals"
## 10 "$ind$contrib" "contributions of the individuals"
## 11 "$call" "summary statistics"
## 12 "$call$centre" "mean of the variables"
## 13 "$call$ecart.type" "standard error of the variables"
## 14 "$call$row.w" "weights for the individuals"
## 15 "$call$col.w" "weights for the variables"pcaB$x[,1]#is what we use as explanatory variable. ## [1] -0.9270989 -0.7787576 0.5811682 0.7230741 -0.1059241 -0.1059241
## [7] -0.1059241 0.5811682 -0.9605731 -0.2877396 -0.7787576 0.7230741
## [13] 0.5811682 -0.9605731 -0.9270989 -0.1059241 -0.7787576 -0.1059241
## [19] 0.3993527 -0.7787576 -0.1059241 0.5811682 -0.9605731 -0.9270989
## [25] 0.5811682 0.5811682 0.3993527 0.3993527 -0.1059241 0.5811682
## [31] -0.0916653 -0.7787576 0.5811682 0.5747328 -0.4218221 0.5811682
## [37] 0.5811682 0.5811682 -0.9605731 0.5811682 -0.7787576 -0.2400066
## [43] 0.5811682 0.7230741 -0.9605731 0.5811682 0.3993527 0.5811682
## [49] 0.5811682 0.3993527#with PC1 as explanatory variable:
Bbsp<-glm(data$Bbuf_pres ~ pcaB$x[,1], family = "binomial")
testZeroInflation(simulateResiduals(Bbsp))
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 0.99871, p-value = 1
## alternative hypothesis: two.sidedsimulateResiduals(Bbsp, plot=T)
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.2562227 0.2492455 0.7444733 0.9873053 0.03247239 0.0819542 0.4597335 0.06559149 0.0448016 0.1499663 0.4141015 0.2128683 0.5981021 0.7153422 0.9464596 0.4942223 0.5403071 0.1554513 0.9880159 0.2774713 ...summary(Bbsp)##
## Call:
## glm(formula = data$Bbuf_pres ~ pcaB$x[, 1], family = "binomial")
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.4910 0.2919 -1.682 0.0925 .
## pcaB$x[, 1] 0.1775 0.4723 0.376 0.7070
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 66.406 on 49 degrees of freedom
## Residual deviance: 66.264 on 48 degrees of freedom
## AIC: 70.264
##
## Number of Fisher Scoring iterations: 4Bb<-glm(data$Bbuf_pres~ data$X1_submerg + data$X1_emerg + data$X1_emerg_Land + data$mean.mud + data$forest.dist + data$MeanTemp + pcaB$x[,1], family = "binomial")
testZeroInflation(simulateResiduals(Bb))
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 0.99512, p-value = 1
## alternative hypothesis: two.sidedsimulateResiduals(Bb, plot=T)
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.2082843 0.2280974 0.7469069 0.9897005 0.03628088 0.06526891 0.5528821 0.06559149 0.04862612 0.1682549 0.426804 0.1978776 0.5380052 0.7645423 0.9301942 0.5498223 0.5339506 0.09838691 0.9897279 0.2477422 ...drop1(Bb)# PCO1 --> other sp. does not explain the probability of presence of B. bufo. ## Single term deletions
##
## Model:
## data$Bbuf_pres ~ data$X1_submerg + data$X1_emerg + data$X1_emerg_Land +
## data$mean.mud + data$forest.dist + data$MeanTemp + pcaB$x[,
## 1]
## Df Deviance AIC
## <none> 60.484 76.484
## data$X1_submerg 1 60.762 74.762
## data$X1_emerg 1 61.076 75.076
## data$X1_emerg_Land 1 61.347 75.347
## data$mean.mud 1 63.448 77.448
## data$forest.dist 1 60.637 74.637
## data$MeanTemp 1 60.518 74.518
## pcaB$x[, 1] 1 60.491 74.491Bb<-glm(Bbuf_pres~ X1_submerg + X1_emerg + X1_emerg_Land +
mean.mud + forest.dist + MeanTemp, data, family = "binomial")
testZeroInflation(simulateResiduals(Bb))
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 0.99538, p-value = 1
## alternative hypothesis: two.sidedsimulateResiduals(Bb, plot=T)
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.2082843 0.2235656 0.7493404 0.9898203 0.03628088 0.06526891 0.5588917 0.06643783 0.04835294 0.1682549 0.421723 0.1978776 0.5530294 0.761028 0.9261278 0.5498223 0.5307723 0.09346756 0.9898502 0.2427874 ...drop1(Bb)# Mean_Temp## Single term deletions
##
## Model:
## Bbuf_pres ~ X1_submerg + X1_emerg + X1_emerg_Land + mean.mud +
## forest.dist + MeanTemp
## Df Deviance AIC
## <none> 60.491 74.491
## X1_submerg 1 60.777 72.777
## X1_emerg 1 61.134 73.134
## X1_emerg_Land 1 61.355 73.355
## mean.mud 1 63.453 75.453
## forest.dist 1 60.656 72.656
## MeanTemp 1 60.524 72.524Bb1<-glm(Bbuf_pres~ X1_submerg + X1_emerg + X1_emerg_Land +
mean.mud + forest.dist,
data, family = "binomial")
simulateResiduals(Bb1, plot=T)
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.2099373 0.2280974 0.7493404 0.9898203 0.03628088 0.06183371 0.5679061 0.06516832 0.04862612 0.166426 0.426804 0.1918813 0.5680537 0.7469708 0.9261278 0.5498223 0.5244157 0.09838691 0.9902171 0.2345293 ...drop1(Bb1)# forest.dist ## Single term deletions
##
## Model:
## Bbuf_pres ~ X1_submerg + X1_emerg + X1_emerg_Land + mean.mud +
## forest.dist
## Df Deviance AIC
## <none> 60.524 72.524
## X1_submerg 1 60.862 70.862
## X1_emerg 1 61.275 71.275
## X1_emerg_Land 1 61.365 71.365
## mean.mud 1 63.465 73.465
## forest.dist 1 60.662 70.662Bb2<-glm(Bbuf_pres~ X1_submerg + X1_emerg + X1_emerg_Land +
mean.mud,
data, family = "binomial")
simulateResiduals(Bb2, plot=T)
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.2082843 0.2160127 0.7396061 0.9895808 0.03668177 0.06281519 0.5498773 0.06559149 0.04780658 0.1636828 0.4395065 0.1813878 0.5530294 0.7469708 0.9234169 0.5405556 0.5403071 0.09838691 0.9892388 0.2295745 ...drop1(Bb2)#submerg## Single term deletions
##
## Model:
## Bbuf_pres ~ X1_submerg + X1_emerg + X1_emerg_Land + mean.mud
## Df Deviance AIC
## <none> 60.662 70.662
## X1_submerg 1 60.917 68.917
## X1_emerg 1 61.305 69.305
## X1_emerg_Land 1 61.371 69.371
## mean.mud 1 63.558 71.558Bb3<-glm(Bbuf_pres~ X1_emerg + X1_emerg_Land +
mean.mud,
data, family = "binomial")
simulateResiduals(Bb3, plot=T)
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.1950599 0.2160127 0.7396061 0.9898203 0.03648132 0.06772263 0.5378582 0.07024637 0.04726022 0.166426 0.4344255 0.1933804 0.5718097 0.7434565 0.9234169 0.5467334 0.5371289 0.1062579 0.9894834 0.2345293 ...drop1(Bb3)#emerg ## Single term deletions
##
## Model:
## Bbuf_pres ~ X1_emerg + X1_emerg_Land + mean.mud
## Df Deviance AIC
## <none> 60.917 68.917
## X1_emerg 1 61.307 67.307
## X1_emerg_Land 1 61.595 67.595
## mean.mud 1 63.801 69.801summary(Bb3) # best model since without emerg and emerg_Land Quant dev is observed.##
## Call:
## glm(formula = Bbuf_pres ~ X1_emerg + X1_emerg_Land + mean.mud,
## family = "binomial", data = data)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.5181 0.3099 -1.672 0.0946 .
## X1_emerg -0.2639 0.4244 -0.622 0.5340
## X1_emerg_Land 0.3119 0.3789 0.823 0.4105
## mean.mud 0.6599 0.4251 1.552 0.1206
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 66.406 on 49 degrees of freedom
## Residual deviance: 60.917 on 46 degrees of freedom
## AIC: 68.917
##
## Number of Fisher Scoring iterations: 4tab_model(Bb3)| Bbuf_pres | |||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| (Intercept) | 0.60 | 0.32 – 1.08 | 0.095 |
| X1 emerg | 0.77 | 0.32 – 1.76 | 0.534 |
| X1 emerg Land | 1.37 | 0.65 – 2.93 | 0.410 |
| mean mud | 1.93 | 0.91 – 4.99 | 0.121 |
| Observations | 50 | ||
| R2 Tjur | 0.105 | ||
# None of the environemental variable affect probability of presence of Bufo bufo. We need to extract the results from comparison between models:
# Extract AIC values
aic_values <- c(AIC(Bb), AIC(Bb1), AIC(Bb2), AIC(Bb3))
names(aic_values) <- c("Bb", "Bb1", "Bb2", "Bb3")
# Calculate ΔAIC
aic_min <- min(aic_values)
delta_aic <- aic_values - aic_min
# Calculate Akaike weights
aic_weights <- exp(-0.5 * delta_aic) / sum(exp(-0.5 * delta_aic))
# Create a summary table
model_comparison <- data.frame(
Model = names(aic_values),
AIC = aic_values,
Delta_AIC = delta_aic,
Akaike_Weight = aic_weights
)
print(model_comparison)## Model AIC Delta_AIC Akaike_Weight
## Bb Bb 74.49068 5.573976 0.03746865
## Bb1 Bb1 72.52442 3.607712 0.10014676
## Bb2 Bb2 70.66152 1.744810 0.25419140
## Bb3 Bb3 68.91671 0.000000 0.60819320Now, we need to represent it graphically:
summ(Bb3)| Observations | 50 |
| Dependent variable | Bbuf_pres |
| Type | Generalized linear model |
| Family | binomial |
| Link | logit |
| χ²(3) | 5.49 |
| Pseudo-R² (Cragg-Uhler) | 0.14 |
| Pseudo-R² (McFadden) | 0.08 |
| AIC | 68.92 |
| BIC | 76.56 |
| Est. | S.E. | z val. | p | |
|---|---|---|---|---|
| (Intercept) | -0.52 | 0.31 | -1.67 | 0.09 |
| X1_emerg | -0.26 | 0.42 | -0.62 | 0.53 |
| X1_emerg_Land | 0.31 | 0.38 | 0.82 | 0.41 |
| mean.mud | 0.66 | 0.43 | 1.55 | 0.12 |
| Standard errors: MLE |
plot_summs(Bb3, colors = "#9E6969", point.size = 1, inner_ci_level = .95)## Le chargement a nécessité le package : broom.mixed
to be able to plot all the species on the same graph, we need to extract the coeff and se of each model. using the following function.
sumbb<-summary(Bb3)
coef<-sumbb$coefficients[,1:2]
dfcoef<-as.data.frame(coef)
dfcoef$variable<-row.names(dfcoef)
dfcoef$species_name<-"Bufo bufo"
coefBB<-dfcoefHyla arborea
length(data$points[data$Harbo_pres==1])## [1] 4PCA:
pcaH<-prcomp(data[,c(11,13:16)])
pcaH## Standard deviations (1, .., p=5):
## [1] 0.6238767 0.5337769 0.4227561 0.3253126 0.2366494
##
## Rotation (n x k) = (5 x 5):
## PC1 PC2 PC3 PC4 PC5
## Bbuf_pres 0.09524347 0.8570632 0.16128109 -0.44887571 0.1699129
## Lvulga_pres -0.16039910 0.2781172 -0.91825270 0.08550326 -0.2154626
## Pridi_pres -0.13833213 0.3171985 0.33981549 0.40499699 -0.7750820
## Pberg_pres -0.67017677 0.2094372 0.10692662 0.48811948 0.5072526
## Rtemp_pres 0.70493130 0.2088413 -0.06239011 0.62363178 0.2581628str(pcaH)## List of 5
## $ sdev : num [1:5] 0.624 0.534 0.423 0.325 0.237
## $ rotation: num [1:5, 1:5] 0.0952 -0.1604 -0.1383 -0.6702 0.7049 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:5] "Bbuf_pres" "Lvulga_pres" "Pridi_pres" "Pberg_pres" ...
## .. ..$ : chr [1:5] "PC1" "PC2" "PC3" "PC4" ...
## $ center : Named num [1:5] 0.38 0.24 0.12 0.34 0.56
## ..- attr(*, "names")= chr [1:5] "Bbuf_pres" "Lvulga_pres" "Pridi_pres" "Pberg_pres" ...
## $ scale : logi FALSE
## $ x : num [1:50, 1:5] -0.957 -0.818 0.652 0.652 -0.148 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : NULL
## .. ..$ : chr [1:5] "PC1" "PC2" "PC3" "PC4" ...
## - attr(*, "class")= chr "prcomp"PCA(data[,c(11, 13:16)])
## **Results for the Principal Component Analysis (PCA)**
## The analysis was performed on 50 individuals, described by 5 variables
## *The results are available in the following objects:
##
## name description
## 1 "$eig" "eigenvalues"
## 2 "$var" "results for the variables"
## 3 "$var$coord" "coord. for the variables"
## 4 "$var$cor" "correlations variables - dimensions"
## 5 "$var$cos2" "cos2 for the variables"
## 6 "$var$contrib" "contributions of the variables"
## 7 "$ind" "results for the individuals"
## 8 "$ind$coord" "coord. for the individuals"
## 9 "$ind$cos2" "cos2 for the individuals"
## 10 "$ind$contrib" "contributions of the individuals"
## 11 "$call" "summary statistics"
## 12 "$call$centre" "mean of the variables"
## 13 "$call$ecart.type" "standard error of the variables"
## 14 "$call$row.w" "weights for the individuals"
## 15 "$call$col.w" "weights for the variables"pcaH$x[,1]#is what we use as explanatory variable. ## [1] -0.95650720 -0.81817507 0.65217646 0.65217646 -0.14799830 -0.14799830
## [7] -0.14799830 0.55693299 -0.97857417 -0.30839740 -0.81817507 0.55693299
## [13] 0.65217646 -0.88333070 -0.86126373 -0.14799830 -0.81817507 -0.14799830
## [19] 0.49177736 -0.81817507 -0.05275483 0.55693299 -0.97857417 -0.86126373
## [25] 0.65217646 0.55693299 0.39653389 0.39653389 -0.05275483 0.55693299
## [31] -0.01800030 -0.81817507 0.55693299 0.51384433 -0.31673154 0.55693299
## [37] 0.55693299 0.55693299 -0.97857417 0.55693299 -0.81817507 -0.15633244
## [43] 0.55693299 0.65217646 -0.88333070 0.55693299 0.49177736 0.65217646
## [49] 0.55693299 0.49177736Hasp<-glm(data$Harbo_pres ~ pcaH$x[,1], family = "binomial")
simulateResiduals(Hasp, plot=T)
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.4132625 0.3776447 0.2521923 0.9902993 0.05011171 0.1226859 0.7511986 0.09479028 0.06829512 0.2286072 0.635125 0.9374768 0.04269034 0.1214264 0.8305683 0.7722223 0.7945693 0.2459673 0.907385 0.4129037 ...summary(Hasp)##
## Call:
## glm(formula = data$Harbo_pres ~ pcaH$x[, 1], family = "binomial")
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -7.807 5.461 -1.430 0.153
## pcaH$x[, 1] 10.723 9.253 1.159 0.247
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 27.877 on 49 degrees of freedom
## Residual deviance: 20.537 on 48 degrees of freedom
## AIC: 24.537
##
## Number of Fisher Scoring iterations: 9Ha<-glm(data$Harbo_pres~data$X1_submerg +
data$X1_emerg +
data$mean.mud +
data$forest.dist + I(data$forest.dist^2)+
data$MeanTemp + pcaH$x[,1], family = "binomial")## Warning: glm.fit: l'algorithme n'a pas convergé## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1simulateResiduals(Ha, plot=T)
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.4132625 0.3776447 0.391603 0.9700596 0.05011171 0.1226859 0.7511986 0.1057927 0.06829512 0.2286072 0.635125 0.3747682 0.0609862 0.1214264 0.8305683 0.7722223 0.7945693 0.2459673 0.9694284 0.4129037 ...testZeroInflation(simulateResiduals(Ha))
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1, p-value = 1
## alternative hypothesis: two.sidedsummary(Ha)##
## Call:
## glm(formula = data$Harbo_pres ~ data$X1_submerg + data$X1_emerg +
## data$mean.mud + data$forest.dist + I(data$forest.dist^2) +
## data$MeanTemp + pcaH$x[, 1], family = "binomial")
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -379.84 94203.24 -0.004 0.997
## data$X1_submerg -12.69 23524.76 -0.001 1.000
## data$X1_emerg -55.34 27075.31 -0.002 0.998
## data$mean.mud -34.12 22540.73 -0.002 0.999
## data$forest.dist -32.83 15099.33 -0.002 0.998
## I(data$forest.dist^2) 29.69 12727.97 0.002 0.998
## data$MeanTemp 23.90 16517.19 0.001 0.999
## pcaH$x[, 1] 527.45 145959.86 0.004 0.997
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 2.7877e+01 on 49 degrees of freedom
## Residual deviance: 4.0644e-08 on 42 degrees of freedom
## AIC: 16
##
## Number of Fisher Scoring iterations: 25drop1(Ha)# l'algorithme ne converge pas --> no matter the variable I try to remove. Cannot say anything about how the other sp. affect the proba of presence of H. arborea. ## Warning: glm.fit: l'algorithme n'a pas convergé
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1## Warning: glm.fit: l'algorithme n'a pas convergé## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1## Warning: glm.fit: l'algorithme n'a pas convergé## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1## Warning: glm.fit: l'algorithme n'a pas convergé## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1## Warning: glm.fit: l'algorithme n'a pas convergé## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1## Single term deletions
##
## Model:
## data$Harbo_pres ~ data$X1_submerg + data$X1_emerg + data$mean.mud +
## data$forest.dist + I(data$forest.dist^2) + data$MeanTemp +
## pcaH$x[, 1]
## Df Deviance AIC
## <none> 0.0000 16.000
## data$X1_submerg 1 0.0000 14.000
## data$X1_emerg 1 0.0000 14.000
## data$mean.mud 1 0.0000 14.000
## data$forest.dist 1 9.3795 23.379
## I(data$forest.dist^2) 1 0.0000 14.000
## data$MeanTemp 1 0.0000 14.000
## pcaH$x[, 1] 1 22.5916 36.592Ha<-glm(Harbo_pres~X1_submerg +
X1_emerg +
mean.mud +
forest.dist + I(forest.dist^2)+
MeanTemp, data, family = "binomial")
simulateResiduals(Ha, plot=T)
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.4000381 0.3519648 0.3665404 0.9972455 0.04971081 0.1192507 0.7271602 0.1045232 0.06610967 0.2030032 0.635125 0.8149314 0.05952253 0.1185122 0.7907011 0.7228001 0.791391 0.1839835 0.8143198 0.3253681 ...testZeroInflation(simulateResiduals(Ha))
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 0.99896, p-value = 1
## alternative hypothesis: two.sidedsummary(Ha)##
## Call:
## glm(formula = Harbo_pres ~ X1_submerg + X1_emerg + mean.mud +
## forest.dist + I(forest.dist^2) + MeanTemp, family = "binomial",
## data = data)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.95002 1.49142 -2.648 0.00809 **
## X1_submerg 0.08516 0.71956 0.118 0.90579
## X1_emerg -0.84708 1.01625 -0.834 0.40454
## mean.mud -0.31057 0.83913 -0.370 0.71130
## forest.dist -1.11363 0.96179 -1.158 0.24692
## I(forest.dist^2) 0.87463 0.92656 0.944 0.34520
## MeanTemp 0.45575 0.99693 0.457 0.64756
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 27.877 on 49 degrees of freedom
## Residual deviance: 22.592 on 43 degrees of freedom
## AIC: 36.592
##
## Number of Fisher Scoring iterations: 6drop1(Ha)#submerged## Single term deletions
##
## Model:
## Harbo_pres ~ X1_submerg + X1_emerg + mean.mud + forest.dist +
## I(forest.dist^2) + MeanTemp
## Df Deviance AIC
## <none> 22.592 36.592
## X1_submerg 1 22.606 34.606
## X1_emerg 1 23.252 35.252
## mean.mud 1 22.747 34.747
## forest.dist 1 24.165 36.165
## I(forest.dist^2) 1 23.448 35.448
## MeanTemp 1 22.799 34.799Ha1<-glm(Harbo_pres~ X1_emerg +
mean.mud +
forest.dist + I(forest.dist^2)+
MeanTemp, data, family = "binomial")
simulateResiduals(Ha1, plot=T)
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.4016911 0.3519648 0.3696732 0.9973652 0.04971081 0.1192507 0.730165 0.1041 0.06610967 0.2030032 0.635125 0.8099295 0.05952253 0.1185122 0.7907011 0.7228001 0.791391 0.1829996 0.8181976 0.3253681 ...drop1(Ha1)#mean.mud## Single term deletions
##
## Model:
## Harbo_pres ~ X1_emerg + mean.mud + forest.dist + I(forest.dist^2) +
## MeanTemp
## Df Deviance AIC
## <none> 22.606 34.606
## X1_emerg 1 23.354 33.354
## mean.mud 1 22.765 32.765
## forest.dist 1 24.265 34.265
## I(forest.dist^2) 1 23.470 33.470
## MeanTemp 1 22.921 32.921Ha2<-glm(Harbo_pres~ X1_emerg +
forest.dist + I(forest.dist^2)+
MeanTemp, data, family = "binomial")
simulateResiduals(Ha2, plot=T)
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.4000381 0.3489437 0.3681068 0.9976048 0.04971081 0.1192507 0.730165 0.1049464 0.06583649 0.2048321 0.635125 0.7699147 0.05927859 0.1189979 0.7740897 0.7320667 0.791391 0.1829996 0.8414638 0.3187617 ...drop1(Ha2)#temp## Single term deletions
##
## Model:
## Harbo_pres ~ X1_emerg + forest.dist + I(forest.dist^2) + MeanTemp
## Df Deviance AIC
## <none> 22.765 32.765
## X1_emerg 1 23.368 31.368
## forest.dist 1 24.448 32.448
## I(forest.dist^2) 1 23.564 31.564
## MeanTemp 1 23.047 31.047Ha3<-glm(Harbo_pres~ X1_emerg +
forest.dist + I(forest.dist^2),
data, family = "binomial")
simulateResiduals(Ha3, plot=T)
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.4000381 0.3459225 0.364974 0.9976048 0.04971081 0.1192507 0.7271602 0.1049464 0.06556331 0.2066609 0.635125 0.8674508 0.05854676 0.1199693 0.7740897 0.7228001 0.791391 0.1623384 0.8065644 0.3270198 ...drop1(Ha3)# forest## Single term deletions
##
## Model:
## Harbo_pres ~ X1_emerg + forest.dist + I(forest.dist^2)
## Df Deviance AIC
## <none> 23.047 31.047
## X1_emerg 1 27.061 33.061
## forest.dist 1 24.521 30.521
## I(forest.dist^2) 1 23.813 29.813summary(Ha3)# best model since without forest dist Quant dev is observed.##
## Call:
## glm(formula = Harbo_pres ~ X1_emerg + forest.dist + I(forest.dist^2),
## family = "binomial", data = data)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.8411 1.4256 -2.694 0.00705 **
## X1_emerg -1.1193 0.6210 -1.802 0.07148 .
## forest.dist -0.9895 0.8586 -1.153 0.24909
## I(forest.dist^2) 0.8155 0.8940 0.912 0.36167
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 27.877 on 49 degrees of freedom
## Residual deviance: 23.047 on 46 degrees of freedom
## AIC: 31.047
##
## Number of Fisher Scoring iterations: 6#hyla arbo tends to be affected by proportion of emerged vegetation. We need to extract the results from comparison between models:
# Extract AIC values
aic_values <- c(AIC(Ha), AIC(Ha1), AIC(Ha2), AIC(Ha3))
names(aic_values) <- c("Ha", "Ha1", "Ha2", "Ha3")
# Calculate ΔAIC
aic_min <- min(aic_values)
delta_aic <- aic_values - aic_min
# Calculate Akaike weights
aic_weights <- exp(-0.5 * delta_aic) / sum(exp(-0.5 * delta_aic))
# Create a summary table
model_comparison <- data.frame(
Model = names(aic_values),
AIC = aic_values,
Delta_AIC = delta_aic,
Akaike_Weight = aic_weights
)
print(model_comparison)## Model AIC Delta_AIC Akaike_Weight
## Ha Ha 36.59158 5.544871 0.03777518
## Ha1 Ha1 34.60555 3.558839 0.10196894
## Ha2 Ha2 32.76496 1.718247 0.25594550
## Ha3 Ha3 31.04671 0.000000 0.60431037Graphic representaiton
summ(Ha3)| Observations | 50 |
| Dependent variable | Harbo_pres |
| Type | Generalized linear model |
| Family | binomial |
| Link | logit |
| χ²(3) | 4.83 |
| Pseudo-R² (Cragg-Uhler) | 0.22 |
| Pseudo-R² (McFadden) | 0.17 |
| AIC | 31.05 |
| BIC | 38.69 |
| Est. | S.E. | z val. | p | |
|---|---|---|---|---|
| (Intercept) | -3.84 | 1.43 | -2.69 | 0.01 |
| X1_emerg | -1.12 | 0.62 | -1.80 | 0.07 |
| forest.dist | -0.99 | 0.86 | -1.15 | 0.25 |
| I(forest.dist^2) | 0.82 | 0.89 | 0.91 | 0.36 |
| Standard errors: MLE |
plot_summs(Ha3, inner_ci_level = .95, colors = "black", plot.distributions = T)
sumbb<-summary(Ha3)
coef<-sumbb$coefficients[,1:2]
dfcoef<-as.data.frame(coef)
dfcoef$variable<-row.names(dfcoef)
dfcoef$species_name<-"Hyla arborea"
coefHA<-dfcoef
coef_final<- full_join(coefBB, coefHA)## Joining with `by = join_by(Estimate, `Std. Error`, variable, species_name)`Lissotriton vulgaris
length(data$points[data$Lvulga_pres==1])## [1] 12PCA
pcaL<-prcomp(data[,c(11:12,14:16)])
pcaL## Standard deviations (1, .., p=5):
## [1] 0.6255488 0.5287522 0.3263076 0.2665902 0.2350240
##
## Rotation (n x k) = (5 x 5):
## PC1 PC2 PC3 PC4 PC5
## Bbuf_pres 0.1688900 0.8570597 0.4461803382 0.1936511 -0.01863562
## Harbo_pres 0.1505903 0.1244790 0.0001832353 -0.7441496 -0.63880266
## Pridi_pres -0.1177651 0.3855602 -0.4560878722 -0.4930396 0.62158714
## Pberg_pres -0.6470553 0.3001270 -0.4741697480 0.2872603 -0.42882144
## Rtemp_pres 0.7185051 0.1059285 -0.6066877081 0.2883298 -0.14603246str(pcaL)## List of 5
## $ sdev : num [1:5] 0.626 0.529 0.326 0.267 0.235
## $ rotation: num [1:5, 1:5] 0.169 0.151 -0.118 -0.647 0.719 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:5] "Bbuf_pres" "Harbo_pres" "Pridi_pres" "Pberg_pres" ...
## .. ..$ : chr [1:5] "PC1" "PC2" "PC3" "PC4" ...
## $ center : Named num [1:5] 0.38 0.08 0.12 0.34 0.56
## ..- attr(*, "names")= chr [1:5] "Bbuf_pres" "Harbo_pres" "Pridi_pres" "Pberg_pres" ...
## $ scale : logi FALSE
## $ x : num [1:50, 1:5] -1.009 -0.892 0.643 0.794 -0.244 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : NULL
## .. ..$ : chr [1:5] "PC1" "PC2" "PC3" "PC4" ...
## - attr(*, "class")= chr "prcomp"PCA(data[,c(11:12, 14:16)])
## **Results for the Principal Component Analysis (PCA)**
## The analysis was performed on 50 individuals, described by 5 variables
## *The results are available in the following objects:
##
## name description
## 1 "$eig" "eigenvalues"
## 2 "$var" "results for the variables"
## 3 "$var$coord" "coord. for the variables"
## 4 "$var$cor" "correlations variables - dimensions"
## 5 "$var$cos2" "cos2 for the variables"
## 6 "$var$contrib" "contributions of the variables"
## 7 "$ind" "results for the individuals"
## 8 "$ind$coord" "coord. for the individuals"
## 9 "$ind$cos2" "cos2 for the individuals"
## 10 "$ind$contrib" "contributions of the individuals"
## 11 "$call" "summary statistics"
## 12 "$call$centre" "mean of the variables"
## 13 "$call$ecart.type" "standard error of the variables"
## 14 "$call$row.w" "weights for the individuals"
## 15 "$call$col.w" "weights for the variables"pcaL$x[,1]## [1] -1.009278032 -0.891512941 0.642937410 0.793527677 -0.244457678
## [6] -0.244457678 -0.244457678 0.474047454 -0.891512941 -0.244457678
## [11] -0.891512941 0.624637721 0.642937410 -0.722622985 -0.840388077
## [16] -0.244457678 -0.891512941 -0.244457678 0.642937410 -0.891512941
## [21] -0.075567723 0.474047454 -0.891512941 -0.840388077 0.642937410
## [26] 0.474047454 0.474047454 0.474047454 -0.075567723 0.474047454
## [31] -0.004117853 -0.891512941 0.474047454 0.675762585 -0.121882945
## [36] 0.474047454 0.474047454 0.474047454 -0.891512941 0.474047454
## [41] -0.891512941 -0.121882945 0.474047454 0.793527677 -0.722622985
## [46] 0.474047454 0.642937410 0.642937410 0.474047454 0.642937410Lvsp<-glm(data$Lvulga_pres ~ pcaL$x[,1], family = "binomial")
simulateResiduals(Lvsp, plot=T)
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.2578758 0.2583089 0.2991847 0.7605268 0.03928758 0.09471353 0.5528821 0.0871732 0.7167617 0.8210369 0.4293445 0.3133062 0.05025263 0.761028 0.5946869 0.5621778 0.5561985 0.1859513 0.9953531 0.2906842 ...summary(Lvsp)##
## Call:
## glm(formula = data$Lvulga_pres ~ pcaL$x[, 1], family = "binomial")
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.1730 0.3369 -3.481 0.000499 ***
## pcaL$x[, 1] -0.4533 0.5302 -0.855 0.392574
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 55.108 on 49 degrees of freedom
## Residual deviance: 54.376 on 48 degrees of freedom
## AIC: 58.376
##
## Number of Fisher Scoring iterations: 4Lv<-glm(data$Lvulga_pres~data$X1_submerg + data$X1_emerg + data$X1_emerg_Land +
data$mean.mud + data$forest.dist + data$MeanTemp + pcaL$x[,1], family = "binomial")
simulateResiduals(Lv, plot=T)
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.332263 0.2462243 0.3696732 0.8885746 0.04730545 0.1035469 0.5889397 0.09394394 0.6012303 0.7963523 0.579234 0.1394138 0.05830281 0.8137424 0.3056492 0.6486667 0.7596083 0.1938222 0.9953531 0.115613 ...testZeroInflation(simulateResiduals(Lv))
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 0.99435, p-value = 1
## alternative hypothesis: two.sidedsummary(Lv)##
## Call:
## glm(formula = data$Lvulga_pres ~ data$X1_submerg + data$X1_emerg +
## data$X1_emerg_Land + data$mean.mud + data$forest.dist + data$MeanTemp +
## pcaL$x[, 1], family = "binomial")
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.44773 0.41601 -3.480 0.000501 ***
## data$X1_submerg 0.20681 0.49568 0.417 0.676509
## data$X1_emerg 1.28867 0.94783 1.360 0.173958
## data$X1_emerg_Land -0.80793 0.51049 -1.583 0.113501
## data$mean.mud 0.08648 0.47946 0.180 0.856855
## data$forest.dist -0.86469 0.49556 -1.745 0.081009 .
## data$MeanTemp 1.32612 0.78721 1.685 0.092071 .
## pcaL$x[, 1] -0.74118 0.65698 -1.128 0.259257
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 55.108 on 49 degrees of freedom
## Residual deviance: 46.540 on 42 degrees of freedom
## AIC: 62.54
##
## Number of Fisher Scoring iterations: 5drop1(Lv)# mean.mud## Single term deletions
##
## Model:
## data$Lvulga_pres ~ data$X1_submerg + data$X1_emerg + data$X1_emerg_Land +
## data$mean.mud + data$forest.dist + data$MeanTemp + pcaL$x[,
## 1]
## Df Deviance AIC
## <none> 46.540 62.540
## data$X1_submerg 1 46.715 60.715
## data$X1_emerg 1 48.731 62.731
## data$X1_emerg_Land 1 49.464 63.464
## data$mean.mud 1 46.572 60.572
## data$forest.dist 1 50.117 64.117
## data$MeanTemp 1 49.813 63.813
## pcaL$x[, 1] 1 47.870 61.870Lv1<-glm(data$Lvulga_pres~data$X1_submerg + data$X1_emerg + data$X1_emerg_Land +data$forest.dist + data$MeanTemp + pcaL$x[,1], family = "binomial")
simulateResiduals(Lv1, plot=T)
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.3471405 0.2492455 0.3696732 0.8885746 0.04730545 0.1040377 0.5889397 0.09267443 0.6086839 0.7932667 0.579234 0.1469091 0.05854676 0.8102281 0.315616 0.6455778 0.7596083 0.1918545 0.9952308 0.1172647 ...drop1(Lv1)#submerg## Single term deletions
##
## Model:
## data$Lvulga_pres ~ data$X1_submerg + data$X1_emerg + data$X1_emerg_Land +
## data$forest.dist + data$MeanTemp + pcaL$x[, 1]
## Df Deviance AIC
## <none> 46.572 60.572
## data$X1_submerg 1 46.742 58.742
## data$X1_emerg 1 48.866 60.866
## data$X1_emerg_Land 1 49.548 61.548
## data$forest.dist 1 50.119 62.119
## data$MeanTemp 1 49.825 61.825
## pcaL$x[, 1] 1 47.932 59.932Lv2<-glm(data$Lvulga_pres~data$X1_emerg + data$X1_emerg_Land +data$forest.dist + data$MeanTemp + pcaL$x[,1], family = "binomial")
simulateResiduals(Lv2, plot=T)
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.3603649 0.2492455 0.3712397 0.8846944 0.04730545 0.1006025 0.6009588 0.08971223 0.627318 0.7901812 0.579234 0.1379147 0.05854676 0.8102281 0.315616 0.6424889 0.7500734 0.1898867 0.9952308 0.107355 ...drop1(Lv2)# pca## Single term deletions
##
## Model:
## data$Lvulga_pres ~ data$X1_emerg + data$X1_emerg_Land + data$forest.dist +
## data$MeanTemp + pcaL$x[, 1]
## Df Deviance AIC
## <none> 46.742 58.742
## data$X1_emerg 1 48.885 58.885
## data$X1_emerg_Land 1 49.623 59.623
## data$forest.dist 1 50.122 60.122
## data$MeanTemp 1 50.406 60.406
## pcaL$x[, 1] 1 48.143 58.143Lv3<-glm(data$Lvulga_pres~data$X1_emerg + data$X1_emerg_Land +data$forest.dist + data$MeanTemp, family = "binomial")
simulateResiduals(Lv3, plot=T)
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.3752423 0.2854994 0.3524427 0.8148501 0.047105 0.1055099 0.5829301 0.08548052 0.7055813 0.7963523 0.604639 0.09144343 0.05537547 0.845371 0.4318955 0.6424889 0.7754997 0.1987416 0.9919291 0.1701163 ...drop1(Lv3)# best model --> sp. does not affect the probability of presence of L. vulgaris. ## Single term deletions
##
## Model:
## data$Lvulga_pres ~ data$X1_emerg + data$X1_emerg_Land + data$forest.dist +
## data$MeanTemp
## Df Deviance AIC
## <none> 48.143 58.143
## data$X1_emerg 1 50.632 58.632
## data$X1_emerg_Land 1 51.149 59.149
## data$forest.dist 1 52.117 60.117
## data$MeanTemp 1 51.303 59.303summary(Lv3)##
## Call:
## glm(formula = data$Lvulga_pres ~ data$X1_emerg + data$X1_emerg_Land +
## data$forest.dist + data$MeanTemp, family = "binomial")
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.3826 0.3974 -3.479 0.000503 ***
## data$X1_emerg 1.1799 0.8163 1.445 0.148337
## data$X1_emerg_Land -0.8353 0.5218 -1.601 0.109414
## data$forest.dist -0.8669 0.4848 -1.788 0.073716 .
## data$MeanTemp 1.2223 0.7487 1.632 0.102582
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 55.108 on 49 degrees of freedom
## Residual deviance: 48.143 on 45 degrees of freedom
## AIC: 58.143
##
## Number of Fisher Scoring iterations: 5Lv<-glm(Lvulga_pres~X1_submerg + X1_emerg + X1_emerg_Land +
mean.mud + forest.dist + MeanTemp,
data, family = "binomial")
simulateResiduals(Lv, plot=T)
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.3603649 0.2854994 0.3493099 0.8187303 0.04690456 0.1074729 0.5528821 0.0901354 0.649679 0.8086946 0.604639 0.0929425 0.05464364 0.8418568 0.4119619 0.6486667 0.7754997 0.2075964 0.9930297 0.1783744 ...testZeroInflation(simulateResiduals(Lv))
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 0.99643, p-value = 1
## alternative hypothesis: two.sideddrop1(Lv)#mean.mud## Single term deletions
##
## Model:
## Lvulga_pres ~ X1_submerg + X1_emerg + X1_emerg_Land + mean.mud +
## forest.dist + MeanTemp
## Df Deviance AIC
## <none> 47.870 61.870
## X1_submerg 1 48.082 60.082
## X1_emerg 1 50.500 62.500
## X1_emerg_Land 1 50.905 62.905
## mean.mud 1 47.932 59.932
## forest.dist 1 52.092 64.092
## MeanTemp 1 50.630 62.630Lv1<-glm(Lvulga_pres~X1_submerg + X1_emerg + X1_emerg_Land +
forest.dist + MeanTemp,
data, family = "binomial")
simulateResiduals(Lv, plot=T)
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.3603649 0.2854994 0.3493099 0.8187303 0.04690456 0.1074729 0.5528821 0.0901354 0.649679 0.8086946 0.604639 0.0929425 0.05464364 0.8418568 0.4119619 0.6486667 0.7754997 0.2075964 0.9930297 0.1783744 ...drop1(Lv1)#submerg## Single term deletions
##
## Model:
## Lvulga_pres ~ X1_submerg + X1_emerg + X1_emerg_Land + forest.dist +
## MeanTemp
## Df Deviance AIC
## <none> 47.932 59.932
## X1_submerg 1 48.143 58.143
## X1_emerg 1 50.622 60.622
## X1_emerg_Land 1 51.014 61.014
## forest.dist 1 52.095 62.095
## MeanTemp 1 50.656 60.656Lv2<-glm(Lvulga_pres~X1_emerg + X1_emerg_Land +
forest.dist + MeanTemp,
data, family = "binomial")
simulateResiduals(Lv2, plot=T)
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.3752423 0.2854994 0.3524427 0.8148501 0.047105 0.1055099 0.5829301 0.08548052 0.7055813 0.7963523 0.604639 0.09144343 0.05537547 0.845371 0.4318955 0.6424889 0.7754997 0.1987416 0.9919291 0.1701163 ...drop1(Lv2)# best model ## Single term deletions
##
## Model:
## Lvulga_pres ~ X1_emerg + X1_emerg_Land + forest.dist + MeanTemp
## Df Deviance AIC
## <none> 48.143 58.143
## X1_emerg 1 50.632 58.632
## X1_emerg_Land 1 51.149 59.149
## forest.dist 1 52.117 60.117
## MeanTemp 1 51.303 59.303summary(Lv2)##
## Call:
## glm(formula = Lvulga_pres ~ X1_emerg + X1_emerg_Land + forest.dist +
## MeanTemp, family = "binomial", data = data)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.3826 0.3974 -3.479 0.000503 ***
## X1_emerg 1.1799 0.8163 1.445 0.148337
## X1_emerg_Land -0.8353 0.5218 -1.601 0.109414
## forest.dist -0.8669 0.4848 -1.788 0.073716 .
## MeanTemp 1.2223 0.7487 1.632 0.102582
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 55.108 on 49 degrees of freedom
## Residual deviance: 48.143 on 45 degrees of freedom
## AIC: 58.143
##
## Number of Fisher Scoring iterations: 5#L. vulgaris tends to be affect by the dist to the wintering habitats. We need to extract the results from comparison between models:
# Extract AIC values
aic_values <- c(AIC(Lv), AIC(Lv1), AIC(Lv2))
names(aic_values) <- c("Lv", "Lv1", "Lv2")
# Calculate ΔAIC
aic_min <- min(aic_values)
delta_aic <- aic_values - aic_min
# Calculate Akaike weights
aic_weights <- exp(-0.5 * delta_aic) / sum(exp(-0.5 * delta_aic))
# Create a summary table
model_comparison <- data.frame(
Model = names(aic_values),
AIC = aic_values,
Delta_AIC = delta_aic,
Akaike_Weight = aic_weights
)
print(model_comparison)## Model AIC Delta_AIC Akaike_Weight
## Lv Lv 61.86963 3.726885 0.09919995
## Lv1 Lv1 59.93204 1.789299 0.26136830
## Lv2 Lv2 58.14274 0.000000 0.63943175Graphical representation:
summ(Lv2)| Observations | 50 |
| Dependent variable | Lvulga_pres |
| Type | Generalized linear model |
| Family | binomial |
| Link | logit |
| χ²(4) | 6.97 |
| Pseudo-R² (Cragg-Uhler) | 0.19 |
| Pseudo-R² (McFadden) | 0.13 |
| AIC | 58.14 |
| BIC | 67.70 |
| Est. | S.E. | z val. | p | |
|---|---|---|---|---|
| (Intercept) | -1.38 | 0.40 | -3.48 | 0.00 |
| X1_emerg | 1.18 | 0.82 | 1.45 | 0.15 |
| X1_emerg_Land | -0.84 | 0.52 | -1.60 | 0.11 |
| forest.dist | -0.87 | 0.48 | -1.79 | 0.07 |
| MeanTemp | 1.22 | 0.75 | 1.63 | 0.10 |
| Standard errors: MLE |
plot_summs(Lv2, inner_ci_level = .95, colors = "black", plot.distributions = T)
sumbb<-summary(Lv2)
coef<-sumbb$coefficients[,1:2]
dfcoef<-as.data.frame(coef)
dfcoef$variable<-row.names(dfcoef)
dfcoef$species_name<-"Lissotriton vulgaris"
coefLv<-dfcoef
coef_final<- full_join(coef_final, coefLv)## Joining with `by = join_by(Estimate, `Std. Error`, variable, species_name)`Pelophylax ridibundus
length(data$points[data$Pridi_pres==1])## [1] 6PCA
pcaPr<-prcomp(data[,c(11:13,15:16)])
pcaPr## Standard deviations (1, .., p=5):
## [1] 0.6269371 0.5159911 0.4089307 0.3066318 0.2479851
##
## Rotation (n x k) = (5 x 5):
## PC1 PC2 PC3 PC4 PC5
## Bbuf_pres 0.1814447 0.85962133 0.3745276 -0.23949634 0.17464118
## Harbo_pres 0.1633600 0.07269459 0.1857509 -0.02906505 -0.96575403
## Lvulga_pres -0.1606287 0.41296865 -0.8734737 -0.12234836 -0.16040537
## Pberg_pres -0.6342752 0.27393148 0.1621125 0.70036667 -0.07656764
## Rtemp_pres 0.7157407 0.10092124 -0.1897073 0.66054078 0.07229882str(pcaPr)## List of 5
## $ sdev : num [1:5] 0.627 0.516 0.409 0.307 0.248
## $ rotation: num [1:5, 1:5] 0.181 0.163 -0.161 -0.634 0.716 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:5] "Bbuf_pres" "Harbo_pres" "Lvulga_pres" "Pberg_pres" ...
## .. ..$ : chr [1:5] "PC1" "PC2" "PC3" "PC4" ...
## $ center : Named num [1:5] 0.38 0.08 0.24 0.34 0.56
## ..- attr(*, "names")= chr [1:5] "Bbuf_pres" "Harbo_pres" "Lvulga_pres" "Pberg_pres" ...
## $ scale : logi FALSE
## $ x : num [1:50, 1:5] -0.863 -0.863 0.669 0.832 -0.229 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : NULL
## .. ..$ : chr [1:5] "PC1" "PC2" "PC3" "PC4" ...
## - attr(*, "class")= chr "prcomp"PCA(data[,c(11:13, 15:16)])
## **Results for the Principal Component Analysis (PCA)**
## The analysis was performed on 50 individuals, described by 5 variables
## *The results are available in the following objects:
##
## name description
## 1 "$eig" "eigenvalues"
## 2 "$var" "results for the variables"
## 3 "$var$coord" "coord. for the variables"
## 4 "$var$cor" "correlations variables - dimensions"
## 5 "$var$cos2" "cos2 for the variables"
## 6 "$var$contrib" "contributions of the variables"
## 7 "$ind" "results for the individuals"
## 8 "$ind$coord" "coord. for the individuals"
## 9 "$ind$cos2" "cos2 for the individuals"
## 10 "$ind$contrib" "contributions of the individuals"
## 11 "$call" "summary statistics"
## 12 "$call$centre" "mean of the variables"
## 13 "$call$ecart.type" "standard error of the variables"
## 14 "$call$row.w" "weights for the individuals"
## 15 "$call$col.w" "weights for the variables"pcaPr$x[,1]## [1] -0.86290331 -0.86290331 0.66855725 0.83191725 -0.22862814 -0.22862814
## [7] -0.22862814 0.48711258 -1.02353198 -0.38925680 -0.86290331 0.65047258
## [13] 0.66855725 -0.84208731 -0.68145865 -0.22862814 -0.86290331 -0.22862814
## [19] 0.50792859 -0.86290331 -0.04718347 0.48711258 -1.02353198 -0.68145865
## [25] 0.66855725 0.48711258 0.32648392 0.32648392 -0.04718347 0.48711258
## [31] 0.03428208 -0.86290331 0.48711258 0.83191725 -0.12634658 0.48711258
## [37] 0.48711258 0.48711258 -1.02353198 0.48711258 -0.86290331 0.03428208
## [43] 0.48711258 0.83191725 -0.84208731 0.48711258 0.50792859 0.66855725
## [49] 0.48711258 0.50792859Prsp<-glm(data$Pridi_pres~ pcaPr$x[,1], family = "binomial")
simulateResiduals(Prsp, plot=T)
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.859183 0.3021157 0.3602748 0.8846944 0.04550143 0.1133618 0.6400212 0.1011378 0.05545564 0.1956878 0.5081 0.3552802 0.05781492 0.09714111 0.9722132 0.6888223 0.6547251 0.2095641 0.8957518 0.3270198 ...summary(Prsp)##
## Call:
## glm(formula = data$Pridi_pres ~ pcaPr$x[, 1], family = "binomial")
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.0678 0.4641 -4.455 8.38e-06 ***
## pcaPr$x[, 1] -0.7130 0.6982 -1.021 0.307
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 36.692 on 49 degrees of freedom
## Residual deviance: 35.632 on 48 degrees of freedom
## AIC: 39.632
##
## Number of Fisher Scoring iterations: 5Pr<-glm(Pridi_pres~X1_emerg +
X1_emerg_Land +
mean.mud + I(mean.mud^2) +
forest.dist +
MeanTemp + pcaPr$x[,1], data, family = "binomial")
simulateResiduals(Pr, plot=T)
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.624488 0.2583089 0.3712397 0.9467782 0.04770634 0.1109081 0.7211506 0.1053696 0.04807976 0.2194629 0.5817745 0.3417886 0.05561942 0.1131694 0.9390046 0.7567778 0.7500734 0.1889029 0.9422844 0.2807745 ...testZeroInflation(simulateResiduals(Pr))
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 0.99233, p-value = 1
## alternative hypothesis: two.sideddrop1(Pr)#emerged vegetation## Single term deletions
##
## Model:
## Pridi_pres ~ X1_emerg + X1_emerg_Land + mean.mud + I(mean.mud^2) +
## forest.dist + MeanTemp + pcaPr$x[, 1]
## Df Deviance AIC
## <none> 28.549 44.549
## X1_emerg 1 28.550 42.550
## X1_emerg_Land 1 28.812 42.812
## mean.mud 1 33.434 47.434
## I(mean.mud^2) 1 30.405 44.405
## forest.dist 1 28.557 42.557
## MeanTemp 1 28.569 42.569
## pcaPr$x[, 1] 1 29.897 43.897Pr1<-glm(Pridi_pres~X1_emerg_Land +
mean.mud + I(mean.mud^2) +
forest.dist +
MeanTemp + pcaPr$x[,1], data, family = "binomial")
simulateResiduals(Pr1, plot=T)
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.624488 0.2583089 0.3712397 0.9467782 0.04770634 0.1109081 0.7211506 0.1053696 0.04807976 0.2194629 0.584315 0.3417886 0.05561942 0.1131694 0.9383269 0.7567778 0.7500734 0.1889029 0.9422844 0.2807745 ...drop1(Pr1)#forest dist## Single term deletions
##
## Model:
## Pridi_pres ~ X1_emerg_Land + mean.mud + I(mean.mud^2) + forest.dist +
## MeanTemp + pcaPr$x[, 1]
## Df Deviance AIC
## <none> 28.550 42.550
## X1_emerg_Land 1 28.845 40.845
## mean.mud 1 33.649 45.649
## I(mean.mud^2) 1 30.406 42.406
## forest.dist 1 28.558 40.558
## MeanTemp 1 28.581 40.581
## pcaPr$x[, 1] 1 29.928 41.928Pr2<-glm(Pridi_pres~X1_emerg_Land +
mean.mud + I(mean.mud^2) +
MeanTemp + pcaPr$x[,1], data, family = "binomial")
simulateResiduals(Pr2, plot=T)
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.6291819 0.2583089 0.3696732 0.9467782 0.04770634 0.1099266 0.7211506 0.1053696 0.04807976 0.2194629 0.584315 0.3402895 0.05537547 0.1131694 0.9376492 0.7567778 0.7532517 0.1889029 0.9422844 0.2758197 ...drop1(Pr2)#Mean temp## Single term deletions
##
## Model:
## Pridi_pres ~ X1_emerg_Land + mean.mud + I(mean.mud^2) + MeanTemp +
## pcaPr$x[, 1]
## Df Deviance AIC
## <none> 28.558 40.558
## X1_emerg_Land 1 28.847 38.847
## mean.mud 1 34.353 44.353
## I(mean.mud^2) 1 30.742 40.742
## MeanTemp 1 28.582 38.582
## pcaPr$x[, 1] 1 29.933 39.933Pr3<-glm(Pridi_pres~X1_emerg_Land +
mean.mud + I(mean.mud^2) +
pcaPr$x[,1], data, family = "binomial")
simulateResiduals(Pr3, plot=T)
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.6338758 0.2583089 0.3712397 0.9467782 0.04790679 0.1069821 0.7211506 0.1053696 0.04862612 0.2185485 0.579234 0.3402895 0.05610731 0.112198 0.9383269 0.7567778 0.7405386 0.1889029 0.9422844 0.2725165 ...drop1(Pr3)#emerg land## Single term deletions
##
## Model:
## Pridi_pres ~ X1_emerg_Land + mean.mud + I(mean.mud^2) + pcaPr$x[,
## 1]
## Df Deviance AIC
## <none> 28.582 38.582
## X1_emerg_Land 1 29.014 37.014
## mean.mud 1 34.569 42.569
## I(mean.mud^2) 1 30.743 38.743
## pcaPr$x[, 1] 1 30.320 38.320summary(Pr3)##
## Call:
## glm(formula = Pridi_pres ~ X1_emerg_Land + mean.mud + I(mean.mud^2) +
## pcaPr$x[, 1], family = "binomial", data = data)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.2974 0.6249 -3.677 0.000236 ***
## X1_emerg_Land 0.3788 0.5774 0.656 0.511875
## mean.mud 2.0199 0.9811 2.059 0.039510 *
## I(mean.mud^2) -0.3983 0.3000 -1.327 0.184346
## pcaPr$x[, 1] -1.0554 0.8243 -1.280 0.200461
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 36.692 on 49 degrees of freedom
## Residual deviance: 28.582 on 45 degrees of freedom
## AIC: 38.582
##
## Number of Fisher Scoring iterations: 6Pr4<-glm(Pridi_pres~X1_emerg_Land + mean.mud + I(mean.mud^2), data, family = "binomial")
simulateResiduals(Pr4, plot=T)#model not converging. 
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.7066312 0.308158 0.3289465 0.9196165 0.04770634 0.1045284 0.7211506 0.1045232 0.05955334 0.2194629 0.604639 0.2923192 0.04756924 0.1146265 0.9647582 0.7567778 0.7659648 0.1918545 0.9151404 0.3270198 ...drop1(Pr4)## Single term deletions
##
## Model:
## Pridi_pres ~ X1_emerg_Land + mean.mud + I(mean.mud^2)
## Df Deviance AIC
## <none> 30.320 38.320
## X1_emerg_Land 1 31.122 37.122
## mean.mud 1 36.139 42.139
## I(mean.mud^2) 1 32.955 38.955We need to extract the results from comparison between models:
# Extract AIC values
aic_values <- c(AIC(Pr), AIC(Pr1), AIC(Pr2), AIC(Pr3))
names(aic_values) <- c("Pr", "Pr1", "Pr2", "Pr3")
# Calculate ΔAIC
aic_min <- min(aic_values)
delta_aic <- aic_values - aic_min
# Calculate Akaike weights
aic_weights <- exp(-0.5 * delta_aic) / sum(exp(-0.5 * delta_aic))
# Create a summary table
model_comparison <- data.frame(
Model = names(aic_values),
AIC = aic_values,
Delta_AIC = delta_aic,
Akaike_Weight = aic_weights
)
print(model_comparison)## Model AIC Delta_AIC Akaike_Weight
## Pr Pr 44.54933 5.967068 0.03243256
## Pr1 Pr1 42.55035 3.968080 0.08811621
## Pr2 Pr2 40.55754 1.975276 0.23866451
## Pr3 Pr3 38.58227 0.000000 0.64078671summ(Pr3)| Observations | 50 |
| Dependent variable | Pridi_pres |
| Type | Generalized linear model |
| Family | binomial |
| Link | logit |
| χ²(4) | 8.11 |
| Pseudo-R² (Cragg-Uhler) | 0.29 |
| Pseudo-R² (McFadden) | 0.22 |
| AIC | 38.58 |
| BIC | 48.14 |
| Est. | S.E. | z val. | p | |
|---|---|---|---|---|
| (Intercept) | -2.30 | 0.62 | -3.68 | 0.00 |
| X1_emerg_Land | 0.38 | 0.58 | 0.66 | 0.51 |
| mean.mud | 2.02 | 0.98 | 2.06 | 0.04 |
| I(mean.mud^2) | -0.40 | 0.30 | -1.33 | 0.18 |
| pcaPr$x[, 1] | -1.06 | 0.82 | -1.28 | 0.20 |
| Standard errors: MLE |
plot_summs(Pr3, inner_ci_level = .95, colors = "Black", plot.distributions = T)
sumbb<-summary(Pr3)
coef<-sumbb$coefficients[,1:2]
dfcoef<-as.data.frame(coef)
dfcoef$variable<-row.names(dfcoef)
dfcoef$species_name<-"Pelophylax ridibundus"
coefPr<-dfcoef
coef_final<- full_join(coef_final, coefPr)## Joining with `by = join_by(Estimate, `Std. Error`, variable, species_name)`Pelophylax bergeri
length(data$points[data$Pberg_pres==1])## [1] 17PCA
pcaPb<-prcomp(data[,c(11:14,16)])
pcaPb## Standard deviations (1, .., p=5):
## [1] 0.5568887 0.4812606 0.4214406 0.2822512 0.2485994
##
## Rotation (n x k) = (5 x 5):
## PC1 PC2 PC3 PC4 PC5
## Bbuf_pres 0.70078214 0.52160158 0.2592054 0.38696777 0.14108400
## Harbo_pres 0.19546022 -0.09325664 0.1415840 0.02890934 -0.96551371
## Lvulga_pres 0.04402947 0.53328906 -0.8038585 -0.19942801 -0.16644566
## Pridi_pres 0.17650508 0.22390988 0.3688108 -0.88371245 0.04172784
## Rtemp_pres 0.66151876 -0.62024361 -0.3613262 -0.16941362 0.13576896str(pcaPb)## List of 5
## $ sdev : num [1:5] 0.557 0.481 0.421 0.282 0.249
## $ rotation: num [1:5, 1:5] 0.701 0.195 0.044 0.177 0.662 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:5] "Bbuf_pres" "Harbo_pres" "Lvulga_pres" "Pridi_pres" ...
## .. ..$ : chr [1:5] "PC1" "PC2" "PC3" "PC4" ...
## $ center : Named num [1:5] 0.38 0.08 0.24 0.12 0.56
## ..- attr(*, "names")= chr [1:5] "Bbuf_pres" "Harbo_pres" "Lvulga_pres" "Pridi_pres" ...
## $ scale : logi FALSE
## $ x : num [1:50, 1:5] -0.508 -0.684 0.678 0.874 -0.684 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : NULL
## .. ..$ : chr [1:5] "PC1" "PC2" "PC3" "PC4" ...
## - attr(*, "class")= chr "prcomp"PCA(data[,c(11:14,16)])
## **Results for the Principal Component Analysis (PCA)**
## The analysis was performed on 50 individuals, described by 5 variables
## *The results are available in the following objects:
##
## name description
## 1 "$eig" "eigenvalues"
## 2 "$var" "results for the variables"
## 3 "$var$coord" "coord. for the variables"
## 4 "$var$cor" "correlations variables - dimensions"
## 5 "$var$cos2" "cos2 for the variables"
## 6 "$var$contrib" "contributions of the variables"
## 7 "$ind" "results for the individuals"
## 8 "$ind$coord" "coord. for the individuals"
## 9 "$ind$cos2" "cos2 for the individuals"
## 10 "$ind$contrib" "contributions of the individuals"
## 11 "$call" "summary statistics"
## 12 "$call$centre" "mean of the variables"
## 13 "$call$ecart.type" "standard error of the variables"
## 14 "$call$row.w" "weights for the individuals"
## 15 "$call$col.w" "weights for the variables"pcaPb$x[,1]## [1] -0.50762714 -0.68413222 0.67816868 0.87362889 -0.68413222 -0.68413222
## [7] -0.68413222 -0.02261346 -0.64010274 -0.64010274 -0.68413222 0.17284676
## [13] 0.67816868 0.06067939 0.19315500 -0.68413222 -0.68413222 -0.68413222
## [19] 0.72219815 -0.68413222 0.01664992 -0.02261346 -0.64010274 0.19315500
## [25] 0.67816868 -0.02261346 0.02141601 0.02141601 0.01664992 -0.02261346
## [31] 0.67816868 -0.68413222 -0.02261346 1.05013397 0.89870323 -0.02261346
## [37] -0.02261346 -0.02261346 -0.64010274 -0.02261346 -0.68413222 0.85467376
## [43] -0.02261346 0.87362889 0.06067939 -0.02261346 0.72219815 0.67816868
## [49] -0.02261346 0.72219815Pbsp<-glm(Pberg_pres ~ pcaPb$x[,1], data, family = "binomial")
simulateResiduals(Pbsp, plot=T)
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.7042843 0.6838435 0.3023175 0.7954489 0.02245004 0.05741701 0.3846137 0.07363174 0.5304207 0.115218 0.8117245 0.2713321 0.05049658 0.7504851 0.9566255 0.3459556 0.8857805 0.1151127 0.8492193 0.6735745 ...summary(Pbsp)##
## Call:
## glm(formula = Pberg_pres ~ pcaPb$x[, 1], family = "binomial",
## data = data)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.7280 0.3181 -2.289 0.0221 *
## pcaPb$x[, 1] -1.1420 0.6046 -1.889 0.0589 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 64.104 on 49 degrees of freedom
## Residual deviance: 60.150 on 48 degrees of freedom
## AIC: 64.15
##
## Number of Fisher Scoring iterations: 4Pb<-glm(Pberg_pres~X1_submerg + X1_emerg + X1_emerg_Land + mean.mud + forest.dist + MeanTemp + pcaPb$x[,1],
data, family = "binomial")
simulateResiduals(Pb, plot=T)
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.5728551 0.6863329 0.2913526 0.8769339 0.01743887 0.07459304 0.291465 0.0930976 0.3366261 0.147223 0.7854535 0.2953173 0.05074052 0.6837135 0.9471373 0.4973111 0.8701678 0.2115319 0.9384067 0.7956905 ...testZeroInflation(simulateResiduals(Pb))
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 0.99338, p-value = 1
## alternative hypothesis: two.sideddrop1(Pb) # MeanTemp## Single term deletions
##
## Model:
## Pberg_pres ~ X1_submerg + X1_emerg + X1_emerg_Land + mean.mud +
## forest.dist + MeanTemp + pcaPb$x[, 1]
## Df Deviance AIC
## <none> 53.358 69.358
## X1_submerg 1 53.579 67.579
## X1_emerg 1 56.597 70.597
## X1_emerg_Land 1 53.556 67.556
## mean.mud 1 55.692 69.692
## forest.dist 1 53.571 67.571
## MeanTemp 1 53.440 67.440
## pcaPb$x[, 1] 1 56.650 70.650Pb1<-glm(Pberg_pres~X1_submerg + X1_emerg + X1_emerg_Land + mean.mud + forest.dist + pcaPb$x[,1],
data, family = "binomial")
simulateResiduals(Pb1, plot=T)
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.5681612 0.6938012 0.2897862 0.8769339 0.01663709 0.07802825 0.2884602 0.09563662 0.3291725 0.1499663 0.7796155 0.308809 0.04976474 0.7083136 0.9498482 0.4880445 0.8709895 0.2085802 0.934529 0.7956905 ...drop1(Pb1) # emerg Land## Single term deletions
##
## Model:
## Pberg_pres ~ X1_submerg + X1_emerg + X1_emerg_Land + mean.mud +
## forest.dist + pcaPb$x[, 1]
## Df Deviance AIC
## <none> 53.440 67.440
## X1_submerg 1 53.665 65.665
## X1_emerg 1 58.010 70.010
## X1_emerg_Land 1 53.602 65.602
## mean.mud 1 55.791 67.791
## forest.dist 1 53.837 65.837
## pcaPb$x[, 1] 1 56.694 68.694Pb2<-glm(Pberg_pres~X1_submerg + X1_emerg + mean.mud + forest.dist + pcaPb$x[,1],
data, family = "binomial")
simulateResiduals(Pb2, plot=T)
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.5399978 0.6888223 0.2788213 0.8536525 0.01543441 0.06870412 0.2944698 0.09225126 0.3440798 0.1490519 0.7591825 0.3133062 0.04854502 0.6766849 0.9525591 0.4787778 0.8594854 0.2085802 0.9384067 0.8074324 ...drop1(Pb2) # submerged## Single term deletions
##
## Model:
## Pberg_pres ~ X1_submerg + X1_emerg + mean.mud + forest.dist +
## pcaPb$x[, 1]
## Df Deviance AIC
## <none> 53.602 65.602
## X1_submerg 1 53.816 63.816
## X1_emerg 1 58.226 68.226
## mean.mud 1 56.040 66.040
## forest.dist 1 54.178 64.178
## pcaPb$x[, 1] 1 56.846 66.846Pb3<-glm(Pberg_pres~X1_emerg + mean.mud + forest.dist + pcaPb$x[,1],
data, family = "binomial")
simulateResiduals(Pb3, plot=T)
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.575202 0.6863329 0.2819542 0.8575327 0.01463262 0.0662504 0.3094938 0.08928906 0.3440798 0.1453942 0.757723 0.308809 0.04854502 0.6837135 0.9545923 0.4756889 0.8594854 0.2046448 0.9384067 0.7956905 ...drop1(Pb3) # forst.dist## Single term deletions
##
## Model:
## Pberg_pres ~ X1_emerg + mean.mud + forest.dist + pcaPb$x[, 1]
## Df Deviance AIC
## <none> 53.816 63.816
## X1_emerg 1 59.208 67.208
## mean.mud 1 56.360 64.360
## forest.dist 1 54.721 62.721
## pcaPb$x[, 1] 1 57.417 65.417Pb4<-glm(Pberg_pres~X1_emerg + mean.mud + pcaPb$x[,1],
data, family = "binomial")
simulateResiduals(Pb4, plot=T)
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.5540795 0.6788646 0.2474931 0.8032094 0.01603575 0.05987073 0.2644219 0.08759638 0.3366261 0.1316778 0.7708585 0.2788275 0.04317823 0.6485706 0.9423932 0.4417111 0.863594 0.1967738 0.907385 0.7839486 ...drop1(Pb4)# bestmodel## Single term deletions
##
## Model:
## Pberg_pres ~ X1_emerg + mean.mud + pcaPb$x[, 1]
## Df Deviance AIC
## <none> 54.721 62.721
## X1_emerg 1 59.997 65.997
## mean.mud 1 57.065 63.065
## pcaPb$x[, 1] 1 57.693 63.693summary(Pb4)##
## Call:
## glm(formula = Pberg_pres ~ X1_emerg + mean.mud + pcaPb$x[, 1],
## family = "binomial", data = data)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.8403 0.3504 -2.398 0.0165 *
## X1_emerg 1.0235 0.4987 2.052 0.0401 *
## mean.mud 0.6400 0.4165 1.537 0.1244
## pcaPb$x[, 1] -1.1073 0.6659 -1.663 0.0963 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 64.104 on 49 degrees of freedom
## Residual deviance: 54.721 on 46 degrees of freedom
## AIC: 62.721
##
## Number of Fisher Scoring iterations: 4We need to extract the results from comparison between models:
# Extract AIC values
aic_values <- c(AIC(Pb), AIC(Pb1), AIC(Pb2), AIC(Pb3), AIC(Pb4))
names(aic_values) <- c("Pb", "Pb1", "Pb2", "Pb3", "Pb4")
# Calculate ΔAIC
aic_min <- min(aic_values)
delta_aic <- aic_values - aic_min
# Calculate Akaike weights
aic_weights <- exp(-0.5 * delta_aic) / sum(exp(-0.5 * delta_aic))
# Create a summary table
model_comparison <- data.frame(
Model = names(aic_values),
AIC = aic_values,
Delta_AIC = delta_aic,
Akaike_Weight = aic_weights
)
print(model_comparison)## Model AIC Delta_AIC Akaike_Weight
## Pb Pb 69.35803 6.637270 0.01860592
## Pb1 Pb1 67.44012 4.719367 0.04854208
## Pb2 Pb2 65.60165 2.880890 0.12171344
## Pb3 Pb3 63.81621 1.095456 0.29719446
## Pb4 Pb4 62.72076 0.000000 0.51394410summ(Pb4)| Observations | 50 |
| Dependent variable | Pberg_pres |
| Type | Generalized linear model |
| Family | binomial |
| Link | logit |
| χ²(3) | 9.38 |
| Pseudo-R² (Cragg-Uhler) | 0.24 |
| Pseudo-R² (McFadden) | 0.15 |
| AIC | 62.72 |
| BIC | 70.37 |
| Est. | S.E. | z val. | p | |
|---|---|---|---|---|
| (Intercept) | -0.84 | 0.35 | -2.40 | 0.02 |
| X1_emerg | 1.02 | 0.50 | 2.05 | 0.04 |
| mean.mud | 0.64 | 0.42 | 1.54 | 0.12 |
| pcaPb$x[, 1] | -1.11 | 0.67 | -1.66 | 0.10 |
| Standard errors: MLE |
plot_summs(Pb4, inner_ci_level = .9, colors = "black", plot.distributions = T)
sumbb<-summary(Pb4)
coef<-sumbb$coefficients[,1:2]
dfcoef<-as.data.frame(coef)
dfcoef$variable<-row.names(dfcoef)
dfcoef$species_name<-"Pelophylax bergeri"
coefPb<-dfcoef
coef_final<- full_join(coef_final, coefPb)## Joining with `by = join_by(Estimate, `Std. Error`, variable, species_name)`Rana temporaria
length(data$points[data$Rtemp_pres==1])## [1] 28PCA
pcaRt<-prcomp(data[,c(11:15)])
pcaRt## Standard deviations (1, .., p=5):
## [1] 0.5465325 0.4866070 0.4235086 0.2588685 0.2377199
##
## Rotation (n x k) = (5 x 5):
## PC1 PC2 PC3 PC4 PC5
## Bbuf_pres 0.625089286 0.69848866 -0.02523248 0.3427587 -0.05706788
## Harbo_pres -0.003738853 0.26551360 0.11132516 -0.6440599 -0.70871860
## Lvulga_pres 0.345240883 -0.19486754 -0.87020121 -0.2881804 0.05037168
## Pridi_pres 0.355978808 0.03901502 0.36656799 -0.5993301 0.61497015
## Pberg_pres 0.602774607 -0.63412979 0.30878437 0.1595581 -0.33724709str(pcaRt)## List of 5
## $ sdev : num [1:5] 0.547 0.487 0.424 0.259 0.238
## $ rotation: num [1:5, 1:5] 0.62509 -0.00374 0.34524 0.35598 0.60277 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:5] "Bbuf_pres" "Harbo_pres" "Lvulga_pres" "Pridi_pres" ...
## .. ..$ : chr [1:5] "PC1" "PC2" "PC3" "PC4" ...
## $ center : Named num [1:5] 0.38 0.08 0.24 0.12 0.34
## ..- attr(*, "names")= chr [1:5] "Bbuf_pres" "Harbo_pres" "Lvulga_pres" "Pridi_pres" ...
## $ scale : logi FALSE
## $ x : num [1:50, 1:5] 0.391 0.035 0.0573 0.0536 -0.5678 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : NULL
## .. ..$ : chr [1:5] "PC1" "PC2" "PC3" "PC4" ...
## - attr(*, "class")= chr "prcomp"PCA(data[,c(11:15)])
## **Results for the Principal Component Analysis (PCA)**
## The analysis was performed on 50 individuals, described by 5 variables
## *The results are available in the following objects:
##
## name description
## 1 "$eig" "eigenvalues"
## 2 "$var" "results for the variables"
## 3 "$var$coord" "coord. for the variables"
## 4 "$var$cor" "correlations variables - dimensions"
## 5 "$var$cos2" "cos2 for the variables"
## 6 "$var$contrib" "contributions of the variables"
## 7 "$ind" "results for the individuals"
## 8 "$ind$coord" "coord. for the individuals"
## 9 "$ind$cos2" "cos2 for the individuals"
## 10 "$ind$contrib" "contributions of the individuals"
## 11 "$call" "summary statistics"
## 12 "$call$centre" "mean of the variables"
## 13 "$call$ecart.type" "standard error of the variables"
## 14 "$call$row.w" "weights for the individuals"
## 15 "$call$col.w" "weights for the variables"pcaRt$x[,1]## [1] 0.39099996 0.03502115 0.05733583 0.05359698 -0.56775346 -0.56775346
## [7] -0.56775346 -0.56775346 0.38026203 -0.22251257 0.03502115 -0.57149231
## [13] 0.05733583 1.00535132 1.01608925 -0.56775346 0.03502115 -0.56775346
## [19] 0.40257671 0.03502115 0.05733583 -0.56775346 0.38026203 1.01608925
## [25] 0.05733583 -0.56775346 -0.22251257 -0.22251257 0.05733583 -0.56775346
## [31] 0.66011044 0.03502115 -0.56775346 0.40957578 1.36133013 -0.56775346
## [37] -0.56775346 -0.56775346 0.38026203 -0.56775346 0.03502115 1.01608925
## [43] -0.56775346 0.05359698 1.00535132 -0.56775346 0.40257671 0.05733583
## [49] -0.56775346 0.40257671Rtsp<- glm(Rtemp_pres~ pcaRt$x[,1],
data, family = "binomial")
simulateResiduals(Rtsp, plot=T)
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.2033251 0.1676742 0.6665984 0.983952 0.01463262 0.03435206 0.2463931 0.3776317 0.03524028 0.08412745 0.2870765 0.5998516 0.4891765 0.07965571 0.5448528 0.2563778 0.3082929 0.07674179 0.9848365 0.1750712 ...summary(Rtsp)##
## Call:
## glm(formula = Rtemp_pres ~ pcaRt$x[, 1], family = "binomial",
## data = data)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.2482 0.2916 0.851 0.395
## pcaRt$x[, 1] -0.8049 0.5492 -1.466 0.143
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 68.593 on 49 degrees of freedom
## Residual deviance: 66.341 on 48 degrees of freedom
## AIC: 70.341
##
## Number of Fisher Scoring iterations: 4Rt<-glm(Rtemp_pres~X1_submerg + X1_emerg + X1_emerg_Land +
mean.mud + I(mean.mud^2) + forest.dist + MeanTemp + pcaRt$x[,1],
data, family = "binomial")
simulateResiduals(Rt, plot=T)
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.332263 0.2432032 0.9342931 0.992455 0.02245004 0.06821337 0.4206712 0.3812086 0.05518246 0.1453942 0.4090205 0.5073173 0.8272215 0.1083123 0.6544879 0.4293556 0.5339506 0.06591923 0.9941302 0.180026 ...testZeroInflation(simulateResiduals(Rt))
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 0.98956, p-value = 1
## alternative hypothesis: two.sideddrop1(Rt)#submerged ## Single term deletions
##
## Model:
## Rtemp_pres ~ X1_submerg + X1_emerg + X1_emerg_Land + mean.mud +
## I(mean.mud^2) + forest.dist + MeanTemp + pcaRt$x[, 1]
## Df Deviance AIC
## <none> 52.815 70.815
## X1_submerg 1 52.819 68.819
## X1_emerg 1 52.868 68.868
## X1_emerg_Land 1 53.663 69.663
## mean.mud 1 52.892 68.892
## I(mean.mud^2) 1 55.137 71.137
## forest.dist 1 55.497 71.497
## MeanTemp 1 53.199 69.199
## pcaRt$x[, 1] 1 56.994 72.994Rt1<-glm(Rtemp_pres~X1_emerg + X1_emerg_Land +
mean.mud + I(mean.mud^2) + forest.dist + MeanTemp + pcaRt$x[,1],
data, family = "binomial")
simulateResiduals(Rt1, plot=T)
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.3339161 0.2432032 0.9342931 0.9923353 0.02265049 0.06772263 0.4206712 0.3776317 0.05545564 0.1453942 0.4090205 0.5073173 0.8272215 0.1073409 0.6544879 0.4293556 0.5339506 0.06591923 0.9941302 0.180026 ...drop1(Rt1)# emerg## Single term deletions
##
## Model:
## Rtemp_pres ~ X1_emerg + X1_emerg_Land + mean.mud + I(mean.mud^2) +
## forest.dist + MeanTemp + pcaRt$x[, 1]
## Df Deviance AIC
## <none> 52.819 68.819
## X1_emerg 1 52.869 66.869
## X1_emerg_Land 1 53.663 67.663
## mean.mud 1 52.895 66.895
## I(mean.mud^2) 1 55.139 69.139
## forest.dist 1 55.568 69.568
## MeanTemp 1 53.212 67.212
## pcaRt$x[, 1] 1 57.143 71.143Rt2<-glm(Rtemp_pres~X1_emerg_Land +
mean.mud + I(mean.mud^2) + forest.dist + MeanTemp + pcaRt$x[,1],
data, family = "binomial")
simulateResiduals(Rt2, plot=T)
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.3355691 0.2432032 0.9342931 0.9922155 0.02265049 0.06870412 0.4296856 0.3812086 0.05545564 0.147223 0.401399 0.5298257 0.8197094 0.1083123 0.6578101 0.4262667 0.5339506 0.06099988 0.9937634 0.1833293 ...drop1(Rt2)# mean temp## Single term deletions
##
## Model:
## Rtemp_pres ~ X1_emerg_Land + mean.mud + I(mean.mud^2) + forest.dist +
## MeanTemp + pcaRt$x[, 1]
## Df Deviance AIC
## <none> 52.869 66.869
## X1_emerg_Land 1 53.678 65.678
## mean.mud 1 52.919 64.919
## I(mean.mud^2) 1 55.183 67.183
## forest.dist 1 55.964 67.964
## MeanTemp 1 53.359 65.359
## pcaRt$x[, 1] 1 57.205 69.205Rt3<-glm(Rtemp_pres~X1_emerg_Land +
mean.mud + I(mean.mud^2) + forest.dist + pcaRt$x[,1],
data, family = "binomial")
simulateResiduals(Rt3, plot=T)
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.2942429 0.2235656 0.9294259 0.9929341 0.02124736 0.07655602 0.3966328 0.4742061 0.04917249 0.1572818 0.391237 0.5598368 0.793417 0.1092838 0.6411988 0.4262667 0.5434854 0.06001601 0.9941302 0.1849809 ...drop1(Rt3)# emerged Land ## Single term deletions
##
## Model:
## Rtemp_pres ~ X1_emerg_Land + mean.mud + I(mean.mud^2) + forest.dist +
## pcaRt$x[, 1]
## Df Deviance AIC
## <none> 53.359 65.359
## X1_emerg_Land 1 54.550 64.550
## mean.mud 1 53.756 63.756
## I(mean.mud^2) 1 55.441 65.441
## forest.dist 1 57.349 67.349
## pcaRt$x[, 1] 1 58.037 68.037Rt4<-glm(Rtemp_pres~mean.mud + I(mean.mud^2) + forest.dist + pcaRt$x[,1],
data, family = "binomial")
simulateResiduals(Rt4, plot=T)
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.2231617 0.2265868 0.8393832 0.9911377 0.01603575 0.05005586 0.3876185 0.4026695 0.05272383 0.153624 0.2972385 0.5798442 0.6882474 0.1058838 0.6677769 0.4262667 0.3082929 0.06395149 0.9955977 0.2378325 ...drop1(Rt4)# mean mud ## Single term deletions
##
## Model:
## Rtemp_pres ~ mean.mud + I(mean.mud^2) + forest.dist + pcaRt$x[,
## 1]
## Df Deviance AIC
## <none> 54.550 64.550
## mean.mud 1 56.119 64.119
## I(mean.mud^2) 1 56.257 64.257
## forest.dist 1 60.775 68.775
## pcaRt$x[, 1] 1 59.729 67.729Rt5<-glm(Rtemp_pres~forest.dist + pcaRt$x[,1],
data, family = "binomial")
simulateResiduals(Rt5, plot=T)
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.2578758 0.2145022 0.7761099 0.9881436 0.01463262 0.0471114 0.3816089 0.3812086 0.04534796 0.1316778 0.248969 0.6973878 0.6206384 0.08985553 0.6544879 0.3428667 0.2606187 0.1023224 0.9910731 0.2692132 ...drop1(Rt5)## Single term deletions
##
## Model:
## Rtemp_pres ~ forest.dist + pcaRt$x[, 1]
## Df Deviance AIC
## <none> 60.542 66.542
## forest.dist 1 66.341 70.341
## pcaRt$x[, 1] 1 62.849 66.849summary(Rt5)##
## Call:
## glm(formula = Rtemp_pres ~ forest.dist + pcaRt$x[, 1], family = "binomial",
## data = data)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.3043 0.3139 0.969 0.3323
## forest.dist 0.7675 0.3426 2.240 0.0251 *
## pcaRt$x[, 1] -0.8772 0.5937 -1.478 0.1395
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 68.593 on 49 degrees of freedom
## Residual deviance: 60.542 on 47 degrees of freedom
## AIC: 66.542
##
## Number of Fisher Scoring iterations: 4We need to extract the results from comparison between models:
# Extract AIC values
aic_values <- c(AIC(Rt), AIC(Rt1), AIC(Rt2), AIC(Rt3), AIC(Rt4), AIC(Rt5))
names(aic_values) <- c("Rt", "Rt1", "Rt2", "Rt3", "Rt4", "Rt5")
# Calculate ΔAIC
aic_min <- min(aic_values)
delta_aic <- aic_values - aic_min
# Calculate Akaike weights
aic_weights <- exp(-0.5 * delta_aic) / sum(exp(-0.5 * delta_aic))
# Create a summary table
model_comparison <- data.frame(
Model = names(aic_values),
AIC = aic_values,
Delta_AIC = delta_aic,
Akaike_Weight = aic_weights
)
print(model_comparison)## Model AIC Delta_AIC Akaike_Weight
## Rt Rt 70.81472 6.2649552 0.01735929
## Rt1 Rt1 68.81935 4.2695814 0.04707843
## Rt2 Rt2 66.86926 2.3194892 0.12481853
## Rt3 Rt3 65.35867 0.8089063 0.26564276
## Rt4 Rt4 64.54977 0.0000000 0.39806112
## Rt5 Rt5 66.54157 1.9918035 0.14703987summ(Rt5)| Observations | 50 |
| Dependent variable | Rtemp_pres |
| Type | Generalized linear model |
| Family | binomial |
| Link | logit |
| χ²(2) | 8.05 |
| Pseudo-R² (Cragg-Uhler) | 0.20 |
| Pseudo-R² (McFadden) | 0.12 |
| AIC | 66.54 |
| BIC | 72.28 |
| Est. | S.E. | z val. | p | |
|---|---|---|---|---|
| (Intercept) | 0.30 | 0.31 | 0.97 | 0.33 |
| forest.dist | 0.77 | 0.34 | 2.24 | 0.03 |
| pcaRt$x[, 1] | -0.88 | 0.59 | -1.48 | 0.14 |
| Standard errors: MLE |
plot_summs(Rt5, inner_ci_level = .95, colors = "black", plot.distributions = T)
sumbb<-summary(Rt5)
coef<-sumbb$coefficients[,1:2]
dfcoef<-as.data.frame(coef)
dfcoef$variable<-row.names(dfcoef)
dfcoef$species_name<-"Rana temporaria"
coefPr<-dfcoef
coef_final<- full_join(coef_final, coefPr)## Joining with `by = join_by(Estimate, `Std. Error`, variable, species_name)`Figure 4
fig.4 <- ggplot(coef_final,
aes(factor(variable,
levels = c("(Intercept)", "I(forest.dist^2)", "forest.dist", "I(mean.mud^2)", "mean.mud", "MeanTemp", "X1_emerg", "X1_emerg_Land", "pcaPb$x[, 1]", "pcaPr$x[, 1]", "pcaRt$x[, 1]"),
labels = c("Intercept", "Distance to Forest^2", "Distance to Forest", "Mud Depth^2", "Mud Depth", "Mean T°C", "Emerged Vegetation", "Emerged Land", "PCA P. bergeri", "PCA P. ridibundus", "PCA R. temporaria")),
Estimate)) +
geom_hline(yintercept=0, lty=2, lwd=1, colour="black") +
geom_errorbar(aes(ymin = Estimate - `Std. Error`, ymax = Estimate + `Std. Error`, colour = species_name), lwd = 1, width = 0) +
geom_point(size=3, aes(colour = species_name)) +
scale_color_manual(values = palette2) +
coord_flip() +
facet_grid(.~species_name) +
guides(colour = "none") +
labs(x="", y="Estimates") +
theme_minimal() +
theme(strip.text = element_text(face = "italic"))
fig.4